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ABSTRACT 



The problem of coupled roll, sway, and yaw stability analysis of submersible 
vehicles is analyzed, with particular emphasis on nonlinear studies. Previous 
results had indicated that a primary loss of stability is through the devel- 
opment of limit cycles. This loss of stability is due to the coupling of roll 
into sway and yaw and cannot be predicted by considering the uncoupled dy- 
namics. In this study, it is shown that the mechanism of loss of stability is 
through bifurcations to periodic solutions. These are characterized as either 
subcritical or supercritical, depending on the sign of a certain nonlinear coeffi- 
cient. Implications of these results to vehicle performance and operations are 



discussed. 
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I. INTRODUCTION 



A. PROBLEM STATEMENT 

The dynamic response of a submersible vehicle operating at the extremes 
of its operational envelope is becoming increasingly important in order to en- 
hance vehicle operations. Traditionally, dynamic stability of motion is studied 
using eigenvalue analysis where the equations of motion are linearized around 
nominal straight line level flight paths [Arentzen & Mandel (I960)], [Clayton 
& Bishop (1982)], [Feldman (1987)]. Directional stability in the horizontal 
plane is normally studied assuming that coupling between sway/yaw and roll 
does not exist. Relaxing this approximation can lead to an oscillatory loss of 
directional stability [Cunningham (1993)] which cannot be predicted by un- 
coupled sway /yaw motions only. This oscillatory loss of stability can generate 
limit cycles in the system, as was confirmed numerically in previous studies 
[Cunningham (1993)]. In order to gain a better understanding of the mecha- 
nism of this type of loss of stability and the stability properties of the resulting 
limit cycles, it is necessary to perform a systematic nonlinear analysis which 
is precisely the scope of this work. 
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B. OBJECTIVES AND OUTLINE 



In this work we examine the problem of stability of motion with controls 
fixed in the horizontal plane, with particular emphasis on the mechanism of 
loss of stability of straight line motion. Coupling between sway/yaw and roll 
is taken into consideration. This has its origins in both inertial and hydro- 
dynamic coupling. We concentrate on an oscillatory loss of stability case, 
where one pair of complex conjugate of eigenvalues of the linearized system 
matrix crosses the imaginary axis. This loss of stability occurs in the form 
of generic bifurcations to periodic solutions [Guckenheimer & Holmes (1983)]. 
Taylor expansions and center manifold approximations are employed in or- 
der to isolate the main nonlinear terms that influence system response after 
the initial loss of stability [Hassard Sz Wan (1978)]. Integral averaging is 
performed in order to combine the nonlinear terms into a design stability co- 
efficient [Chow Sz Mallet-Paret (1977)]. Special attention is paid to the study 
of the quadratic drag terms as they constitute some of the main nonlinear 
terms of the equations of motion. The difficulty associated with the nons- 
moothness of the absolute value nonlinearities is dealt with by employing the 
concept of generalized gradient [Clarke (1983)], a technique which was utilized 
in [Papadimitriou (1994)]. This has the advantage of keeping the linear terms 
constant, unlike the linear/cubic approximation typically used in ship roll mo- 
tion studies [Dalzell (1978)], where the linear damping coefficient is a function 
of the assumed amplitude of motion. 

Vehicle modeling in this work follows standard notation [Gertler Sz Hagen 
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(1976)], [Smith et al (1978)], and numerical results are presented for a variant 
of the Swimmer Delivery Model used in [Cunningham (1993)] for which a set 
of hydrodynamic coefficients and geometric properties is available. Although 
the main results and conclusions of this work are derived for a submerged 
vehicle, similar techniques can be applied to surface ships as well. 
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II. EQUATIONS OF MOTION 



A. COORDINATE SYSTEM 

In our analysis we are going to use a moving coordinate system ( x,y,z ), 
attached on the vehicle. The origin of this system coincides with the center 
of buoyancy, B. The a: -axis is attached to the longitudinal plane of symmetry 
for the vehicle, the y-axis is positive starboard, and the 2 -axis is positive 
downwards. All main symbols used in the development of the equations of 
motion in this work are summarized in Table 1. 



B. GENERAL FORM OF THE EQUATIONS OF MOTION 

The equations of motion for a submersible vehicle in the horizontal plane 
are written as follows: 



Sway equation: 



m[v + Ur — wp + xc(pq + r) — yc(p 2 + r 2 ) + za{qr — p)] = 
Ypp + Yj-r + Y pq pq + Ygr-qr + Y V U v + Y vw vw + 

Ys r U 2 6 r + Yi,v + YpUp + Y r Ur + Y vg vq + Y wp wp + Y wr wr + 
(W — B ) cos 6 sin <j> — 

r x n ose r 0 , ol (v + xr) 

/ \C Dy h{x){y + xrf + C Dz b(x){w - xq) 2 — - - dx . 

x ta.il ^ U c f{X) 



( 1 ) 
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Yaw equation: 



IzzT + ( Iyy ~ Ixx)pq ~ Ixy(p 2 ~ <? 2 ) ~ Iyz(pr + <?) + ^xz(? r - p) + 

m[xc(v + Ur — u>p) — Pg(^ - vr + wg)] = 

Npp + Nj-r + Npgpq + N qr qr + N V U v + N vw vw + 

N( r U“5 r + NyV + NpU p + N r U r + N vq vq + N wp wp + N wr wr + 
(xqW - x B B)cosd sirup + (vgW ~ y B B) sin0 + U 2 N pIop - 

fZnose r n oT (v XT 1 

/ + xr) 2 + C Dz b(x)(w - xq) 2 - - J xdx . 

^tail 1 J U c f\X) 



(2) 



Roll equation: 

IxxP + (Izz - Iyy)qr + I xy (pr - q) - / yz (tf 2 - r 2 ) - / zz (pg 4- r) + 
m [PG(^ — Uq + vp) — 2 g(^ + Ur — wp )] = 

B pP "t~ B j-T “1“ B p q pq H“ i K uU V K-uujV'W 

KpropU 2 + Ki ,v + KpUp + K r U t + K vq vq + K wp wp + K wr wv + 

(vgW — y B B) cosO cos <p — (z G W — z b B) cos 6 sirup ■ ( 3 ) 



The rotational velocity equation around the x-axis is simply, 

<P = P , 

while U c f denotes the cross-flow velocity, 

U c f(x ) = + xr) 2 + (re — x</) 2 . 



( 4 ) 
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b(x) 


local beam of the hull 


C D 


quadratic drag coefficient 


6r 


rudder deflection 


h(x) 


local height of hull 


lyyi Izz) 


vehicle mass moments of inertia about body axes 


(/ X y, Iyz , I Z x ) 


cross products of inertia 


(■ K,M,N ) 


moment components along the three axes 


m 


vehicle mass 


(p,q,r) 


rotational velocity components along the body axes 




Euler angles 


u 


constant vehicle speed along the x-axis 


(u,V,w) 


translational velocities about ( x,y,z ) axes 


(x,y,z) 


distances along the three body axes 


C X,Y,Z ) 


force components along the body axes 


(x G , y G , zg) 


coordinates of the center of gravity 


Vb, zb) 


coordinates of the center of buoyancy 


Znose 


fore coordinate of vehicle body 


^tail 


aft coordinate of vehicle body 



Table T. Nomenclature 



C. SIMPLIFICATIONS 

Before we proceed with the analysis, we must simplify the above equations 
of motion in order to reflect the fact that we are analyzing motion in the 
horizontal plane only. The simplifications that we employ are the following: 



• Velocity, w, and acceleration, w, in the 2 -direction are zero. 



• Acceleration in the longitudinal direction, u, is zero. 

• Rotational velocity, q , and acceleration, q, in the y-direction are zero. 

• Pitch angle, 9 , is zero. 



7 




D. SIMPLIFIED EQUATIONS OF MOTION 



After the application of the above simplifying assumptions, the equations 
of motion become: 

Sway equation: 



(m - Yii)v + {mxG - Yj.)r - ( mza + Y p )p = 



Y v Uv + (Y r U - mU)r 4- Y p Up + yGP 2 + ycr 2 + (W - B) sincf) + Ys T U 2 8 r - 



[ Xnose r w-w- . --n2 {v±xr) 

/ C DyhyX)(y + XT) dx . 

d Xtail U 



U c f(x) 



(5) 



Yaw equation: 

{mxG - Nij)v + ( I zz - Nr)r - ( I xz + N p )p = 

N V U v + ( N r U — mxQU)r + N p Up + I xy p 2 + I yz pr + ycvr + 

(x G W - x b B) sin <f> + N 6r U 2 6 r + U 2 N pTOp - 

Panose (y -L. xr} 

/ C Dy h(x)(v + xr) 2 ^———^ xdx . (6) 

v/ltail UcfyX) 

Roll equation: 

( -mz G - Ky)v + ( I xz - Kf)r + ( I xx - K p )p - 

K V U v + ( K t U + mzGU)r + K p Up — I xy pr — I yz r 2 — mycvp + 

(ycW - y B B) cos $ - (zqW - z b B) sin<£ 4- U 2 K pTOp . (7) 

Roll rate: 

= P ■ (8) 
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III. LINEAR ANALYSIS 



A. LINEARIZATION 

The simplified equations of motion can be written in matrix form as follows: 

Ax = Bx 4- g(x) , (9) 



where the state vector x is defined as, 



x = 



V 

r 

P 

4> 
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and the state matrices are, 



and 



A = 



m — Yy mxc — Yr 
mxo — N{, I zz — Nr 
mzQ Ny Ixz Kf 
0 0 



-mzG -Yp 0 
-Ixz ~Np 0 
0 
1 



I XT Nn 



' Y V U Y r U - mU Y P U 0 
N V U —mx G U + N r U N P U 0 
K V U mz G U + K r U KpU 0 
0 0 10 



The term g(x) contains all the nonlinear terms, 



< 7 i = VGP 2 + yar 2 + (W - B) sin <j> + Yg r U 2 6 r - 



CdJ 

J : 


rX nose 

h(x){y + xr)\v + xr\ dx , 

^tail 


(10) 


92 = IxyP 2 + lyzpr + VGvr -f {xqW — xbB) sm<f> + Ng T U 2 5 r — 

f^nose 


(11) 


C Dy I 

J : 


h(x)(v + xr)\v + xt \ x dx , 

^tail 


93 = ~ IxyPr ■ 


- Iyzl~ 2 - mycvp + U 2 K prop + 
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(y G W - y B B) cos <t> - (z G W - z b B) sin 4> , - 



( 12 ) 



94 



0 . 



(13) 



We want to linearize the nonlinear terms about a nominal point 




After applying Taylor series expansion of the non-linear terms about the nom- 
inal point xo and keeping only the linear components, the linearized equations 
of motion are written in matrix form as: 



performed in the next section in order to assess the dynamic stability of the 
vehicle. 

B. LOSS OF STABILITY 

Stability of the linearized system depends on the location of the four eigen- 
values of the system, in other words the roots of the characteristic equation: 



Ax = B'x , 



(14) 



where 



A' = A , 



and 



' Y V U Y r U -mU Y P U W - B 
N V U -mx G U + N r U N P U x G W - x b B 
K V U mz G U + K r U K P U -z G W + z b B 
0 0 10 



Equation (14) is our linearized system. Eigenvalue analysis for this system is 



det(B' - A A') = 0 



(15) 
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or, in polynomial form, 



AX 4 + BX 3 + CX 2 + DX + E = 0 . (16) 

The coefficients of the polynomial equation (16) are given by: 

A = (Ixx Np)[(m Yy)(I zz N+) (mx G Yf)(mx G — TV^,)] 

+ {piZG + Ki >)[(— -fz 2 + Nr)(mzG + Yp ) + ( Nj, + I xz )(mx G — Yr)) 

“ (Kr)[(-mx G + lVi)(mz G + Yp) + ( N p + I xz )(m - Y*)] , (17) 

B = K V U[(-I ZZ + Nr)(mzG + Yp) + (Np + I xz )(mx G - Yj.)] 

- ( mz G + Ky)[(I zz - N,)(Y P U) - (UN r - Umx G )(mz G + Yp) 

+ (Np + I xz ) (UY r - Um) - ( NpU ) (mx G - Yr)] 

- (I xx - Kp)[(m - Yi,)(UN r - U mx G ) + (Y V U)(I ZZ - N+) 

- (mx G — Yf)(N v U) — (mx G — Ny)(UY r — Um)] 

- K p U[(m - Yy)(I zz - Nj.) - (mx G - Yr)(mx G - N*)] 

+ (. Kr - I xz )[(mx G - Nv)(YpU) - ( N v U)(mz G + Yp) 

+ ( Np + I xz )(Y v U)-(N p U)(m-Yv )] 

- ( mz G U + K r U)[(-mx G + Ny)(mz G + Yp) 

+ (Np + I xz )(m-Yy)] (18) 

C = (I xx - Kp)[(Y v U)(UN r -mUx G ) - (UY r - Um)(N v U)] 

+ K p U[(m - Yi)(UNr - mUx G ) + (Y V U)(I ZZ - N+) 

- (mx G — Yf)(N v U) — (mx G — Ny)(UY r — Um)] 

+ (mz G U + K r U)[(mx G — Ny)(Y p U) — (N v U)(mz G + Yp) 

+ (Np + I XZ )(Y V U) - (N p U)(m - Yy)] 
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-1- ( z G W - z B B)[(m - Yy){I zz - Nr) - (mx G - Yf)(mx G - N o)] 

- ( mz G + Ky)[(x B B - x G W)(mx G - Y+) — (UN r — Umx G )(Y p U ) 

+ (UY r - Um)(N p U)} 

- K V U[(I ZZ - Nr)(Y p U) - ( UN r - Umx G ){mz G + Y p ) 

+ {N p + I xz )(UY r - Um) - (N p U)(mx G - Y f )] 

+ (Kf - I X z)[(xbB - x G W){m - Yy) 

- (N V U)(Y P U) + (N P U)(Y V U)] , (19) 

D = ( mz G U + K r U)[(x B B — x G W){m — Yy) 

- (N v U)(Y p U) + {N P U)(Y V U)\ 

- (Kr - I xz )(x b B - x G W)(Y v U) 

+ (mz G + Ky)(x B B - x G W)(UY r - Um) 

- K V U [{x b B - x G W)(mx G - Yf) 

- (UN r - mUx G )(Y p U) + (UY r - Um)(N p U)\ 

- K p U[(Y v U)(UN r - mUx G ) - ( UY r - Um)(N v U )] 

- (z G W - z B B)[(m - Yy)(UN T - mUx G ) + (Y V U)(I ZZ - N,) 

- (mx G - Yr)(N v U) - (mx G - Ny)(UY r - Um)] , (20) 

E = (z G W - z B B)[(Y v U)(UN r - mUx G ) - (UY r - Um)(N v U)] 

+ ( K v U)(x b B - x G W)(UY r - Um) 

- (mz G U + K t U)(x b B - x G W)(Y v U) (21) 

We can examine the stability of the system by utilizing Routh’s criterion. 
Application of this criterion to the characteristic equation (16), reveals that 
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the following two conditions must be satisfied in order to ensure that all roots 
of (16) have negative real parts: 



BCD - AD 2 - EB 2 > 0 , 


(22) 


E > 0 . 


(23) 



If E is less than zero, one real root of (16) becomes positive and the sytem will 
become unstable in a divergent manner (Guckenheimer and Holmes, 1983). 
This is the case of a directionally unstable ship which is well known in the 
literature (Clayton and Bishop, 1982). If, however, condition (22) is violated, 
the system will exhibit an oscillatory motion due to the presence of complex 
conjugate roots with positive real parts. This form of instability is caused by 
the coupling of roll into sway and yaw and is further analyzed in this work. 

In order to compute the limiting case of loss of stability, we consider equa- 
tion (24), 

BCD - AD 2 - EB 2 = 0 . (24) 

The result of this equation will produce the limiting value of zq as a function 
of x q for loss of stability. A curve of this functional form, 

Z G = f{x G ) , 

will be our locus of loss of stability. After some algebra, we can express the 
coefficients of equation (24) in the following form: 

A = A\Zq + A 2 ZG + A3 , (25) 
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where 



Ai = -m 2 (I zz - Nf) 

A 2 = —mYp(I zz — Nr) — mKjj(I zz — N+) + m(Np + I xz )(rriXG — Yr) 

“I - m{Kf ^xz){n^ , '^g ^ v) 

^4.3 {I xx Np)(rn Yy)(I zz N^) (/n K p)(jnxG Yr){mxQ Ny) 

- KyYp(I zz - Nr) + Ky(Np + I xz )(mx G - Y+) 

+ Yp(K, - I xz )(mx G - Ny) - (K, - I xz )(Np + I xz ){m - Yy) 

B = B\z“g + B 2 z g + , ( 26 ) 

where 

J5j = m 2 (UN r — Unixc) + rn 2 U(mxG — Ny) 

B 2 = -m{K v U)(I zz - Nr)-m{I zz -Nr){Y p U) 

+ mYp(U N T — U mxc) + mKy(U N r — Umxc) 

- m(Np 4 - I xz )(UY r — Um) + m(N p U)(mxG — Yf) 

- m(Kr — I XZ )(N V U) + mUYp(mxG — Ny) 

- mU(Np + I xz )(m —Yy) + mUK r (rnxG — Ny) 

Bz = -Yp{K v U)(I zz -Nr) + {K v U)(Np + I xz )(mx G -Yr) 

- Ky{I zz - Nr)(YpU) + KyYp{UNr-Umx G ) 

- Ky(N p + I xz ){UY r - Um) + Ky{NpU)(mx G - Y+) 

- {I xx ~ K p )(m - Yy){UN r - Umx G ) - (/** - Kp)(Y v U)(I zz - N,) 

+ {Ixx - Kp)(mxG - Yr)(NyU) + ( I xx - Kp)(mxG - Ny)(UY r - Um) 
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- (K p U)(m - Yv)(I zz - Nr) -f ( K p U)(mx G - Yf){mxG - Ny) 

+ (K, - I xz ){mx G - Ny)(Y p U) - Yp(Kf - I XZ )(N V U) 

+ ( N p + I XZ )(Y V U)(K, - I xz ) - (K, - I xz ){N p U){m - Y„) 

+ U K r Yp(mxG - Ni >) 

C = C\Zq + C 2 zg + C 3 > (27) 

where 

C x = —m 2 U (N V U) 

C 2 = mU(mx G - Ny)(Y p U) - mUYp(N v U) - mUK r (N v U ) 

+ mU(Np + I xz )(Y v U)-mU{NpU){m-Yy) 

+ W(m — Yy){I zz — Nr) — W ( miG — Yf)(mxc — Ny) 

- m(X B B — xcW)(mxG — Y+ ) + m{U N r — UmxG)(Y p U) 

- m(UY r - Um)(N p U) +m(K v U)(UN r - Umx G ) 

C 3 = (I xx - Kp)(Y v U)(UN r - Umx G ) - (I xx - Kp)(UY r - Um)(N v U) 

+ ( KpU)(m - Yy){UN r - Umx G ) + ( K P U){Y V U){I ZZ - Nf) 

+ ( K p U)(mx G - Yr)(N v U) - ( K p U)(mx G - Ny)(UY r - Um) 

+ UK r (mx G ~ Ny)(Y p U) - UK r Yp(N v U) 

+ UK r {Np + I xz )(Y v U) - UK r {NpU)(m - Y„) 

- Ky{X B B - x G W){mx G - Yr) + Ky{UN r - Umx G )(Y p U ) 

- Ky(UY r - U m)(N p U) - ( K V U){I ZZ - Nr)(Y p U) 

+ Yp(K v U)(UN r - Umx G ) - ( K V U){N P + I xz ){UY r - Um) 

+ ( K v U)(NpU){mx G - Y,) + {K, - I xz )(X b B - x G W)(m - Y„) 
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- (K, - Ixz)(N v U)(Y p U) + (Kr ~ Ixz)(N p U)(Y v U ) 

D = D\zq + D 2 , (28) 

where 

Di = mU{x s B -x G W)(m -Yj - mU(N v U){Y p U) 

+ mU (N P U)(Y V U) + m(x B B x G W)(UY r - Urn) 

- W (m — Yi)(UN r - Umx G ) ~ W{Y V U){I ZZ - Nf) 

+ W(mx G - Yf)(N v U ) + W(mx G - Ni)(UY r - Um) 

D 2 = U K r (x B B — x G W)(m — Yy) — U K r (N v U)(Y p U) 

+ UK r {N p U)(Y v U) - {Kr - I xz )(x b B - x G W){Y v U) 

+ K^xbB - x G W){UY r - Um) - {K v U){x b B - x G W){mx G - y f ) 

+ {K v U){UN r - Umx G )(Y p U) - ( K v U){UY r - Um){N p U) 

- ( K p U)(Y v U){UN r - Umx G ) + {K P U){UK T - Um){N v U) 

and 

E = E\z g + E 2 , (29) 

where 

Ei - W{Y v U)(UN r -Umx G )-W{UY r -Um){N v U) 

- mU(x B B -x G W){Y v U) 

E 2 = (K v U){x b B -x G W)(UY r -Um)-UK r {x B B -x G W)(Y v U) 

If we apply the stability criterion (24) utilizing expressions (25) through 
(29), we get a fifth order polynomial equation in the metacentric height z G of 
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Figure 1: Critical value of zq versus xq for U = 5 ft/sec 




Figure 2: Critical value of zq versus xq for different values of U (ft/sec) 
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the following form, 



F$ z g + F 4 ZQ + F 3 4 + F^zq + F\Zq + Fq = 0 , (30) 

where, 

Fq = BzCzDi — AzD\ — E2BI , 

Fi = BzCzDi + (.B 3 C 2 + B 2 Cz)D 2 — 2E 2 B 2 Bz ~ F\B^ — 

2D 2 D\Az — D 2 A 2 , 

F 2 = —E 2 {B 2 -+- 2B\Bz) — 2E\B 2 Bz + ( BzCz 4- B 2 Cz)D\ + 

(-B 3 C 1 + B 2 C 2 + B\Cz)D 2 — D\Az — 2D 2 D\A2 — D\A\ , 

Fz = — D\A2 — 2D 2 D\A\ — 2EzB\B 2 — Ei^B^ + 2B\Bz ) + 

Di(BzC\ + B 2 C 2 + B\Cz ) + D 2 (B 2 C\ + BiC 2 ) , 

F4 = D\(^B 2 0\ + BiC 2 ) + B\C\D2 — D^Ai — E2B1 — 2E\B\B2 , 

F 5 = -Erfl + B&Di . 

Using the corresponding Matlab program shown in Appendix A, we can solve 
equation (30) using a typical constant forward velocity U = 5 ft/sec. The value 
of xc is given in the range from 0 to 2 ft. In this way we get five solutions in 
zg for a given value of xq- Investigating these solutions, we can see that only 
one satisfies the second criterion (23). This solution corresponds to the locus 
of loss of stability. For values of zq greater than the critical value, the system 
is stable, and for values less than its critical value, the system is unstable, see 
Figure 1. The calculations are repeated for a range of forward speeds U from 
2 to 8 ft/sec. The results are shown parametrically in Figure 2. We observe 
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an upwards movement of the curves as the speed is increased. This means 
that the system experiences a tendency to become less stable at higher speeds 
which, in turn, calls for higher metacentric heights to ensure stability. 
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IV. NONLINEAR ANALYSIS 



A. INTRODUCTION 

From the linearized analysis on the loss of stability that was done in the 
previous chapter, we can see that as a specific parameter of the system, such as 
( xq , zg), is varied, it is possible to pass from a region of stability to a region of 
instability. This case of loss of stability is associated with one pair of complex 
conjugate eigenvalues of the system crossing the imaginary axis. This loss of 
stability is usually accompanied by self sustained oscillations, and it is called 
Hopf Bifurcation. 

There are two cases of Hopf Bifurcations, supercritical and subcritical. In 
the supercritical case, limit cycles are created just after the loss of stability. 
A limit cycle is a constant amplitude oscillatory motion of our system. This 
amplitude is usually larger as we move away from the bifurcation point. When 
the limit cycles have small amplitudes the situation is not very critical for our 
vehicle and is not much different from stable conditions. In the subcritical 
case we may have convergence to limit cycles, even before our system looses 
its stability. Furthermore, the limit cycle amplitudes are considerably higher. 

The analysis that will follow will be performed in order to verify the exis- 
tence of the limit cycles and to find out which case of Hopf Bifurcations we 
have in our model. We will also examine the stability of the resulting limit 
cycles. This analysis is necessary because we can predict the behavior of a 
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vehicle when some of its parameters change and we have operation in the 
proximity of a bifurcation point. The results of the non linear analysis will 
be verified by a numerical simulation of the system’s motion, both above and 
below the bifurcation point. 



B. THIRD ORDER EXPANSIONS 



From the linearization procedure we see that our system is written in the 
form of matrix equation (14), where we have ignored the non linear terms. 
If we take into account the non linear terms up to third order, equation (14) 
becomes: 

A'x = B'x + g(x) , (31) 

where, 

9i(x) ' 

9i(x) 

9i(x) 

9a{x) . 

Keeping terms up to third order the vector of non linear terms can be written: 



g(x) = 



g(x) = # (2) (x) + # (3) (x) + c(x) , (32) 



where g^ 2 \x) contains the second order non linear terms, g^ z \x) contains the 
third order non linear terms, and c(x) contains the constant terms. The cross 
flow integrals can be written as follows: 



Iv = 

It = 



fZn. 

CdJ 

J x ta.i 
n. 

c dJ 



X n ose 



^tail 

^nose 



h(x)(v + xr)\v + xr| dx 
h(x)(v + xr)\v + xr |x dx 
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( 33 ) 



The second order non linear terms are: 



where 



g (2) (x) 



9 i 2) (z) 
92 \x) 
9z\x) 
9 a\x) 



5 



9\ = UGV 2 + VCr 2 ~ li 2) 

92 } = IxyV 2 + lyzpr + y G vr - I {2) 

9: J 2) = ~ l x y pr - I yz r 2 - my G vr - (y G W - y B B)<f > 2 



The third order non linear terms are: 



where, 




?i 3) (z) 
92 \x) 
93 °(z) 
54 3) (^) 



<4 3) 

4 3) 



„( 3 ) 

94 



-4 3) - \{w - b)4> 3 

-4 3) - kx G w - x b B)4> 3 

0 

Uz G W - z B B)<t> 3 

0 

0 . 



The constant terms are: 



c(x) 



’ ci(x) ' 

C 2 (x) 

C 3 {x) 

. c a{x) _ 
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where, 



ci(x) = Y 6r U 2 6 T 
c 2 {x) = Ns r U 2 6 r + U 2 N wov 

C 3 {x) — U 2 Kprop + VgW — y B B 
C4(x) = 0 . 



In order to get the second and third order non linear terms of the cross flow 
integrals, we must expand in Taylor series a function of the form: 

/(?) = ?l?l (34) 



about a nominal point Co : 



ClC! = ColCol + 2|Co|(C - Co) + sign(Co)(C - Co) 2 + / (3) (C) • (35) 



The sign function in equation (33) can be approximated by: 

sign(Co) = lim tanh(— ) , 

7— >0 7 

where the approximation gets better as 7 gets smaller. 

If we choose Co = 0 as our nominal point, equation (35) becomes: 



(36) 



el«l = 7 -« 3 

67 



In our case C is (v + xr) so we have: 



( 7 ; + xr) lu + xr\ = — (v + xr) 3 
67 



(37) 



(38) 



or 



(?; + xr)\v + xr \ = — ( v 3 + x 3 r 3 + 3 v^xr + 3 x 2 r 2 v) 
67 



( 39 ) 
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Using equation (39) the cross flow integrals become: 

/nr 

I v = -^-{EqV 3 + 3EiV 2 r + 3E 2 vr 2 + E 3 r 3 ) (40) 

07 

I r = -^{E lV 3 + 3E 2 v 2 r + 3E 3 vr 2 + E 4 r 3 ) (41) 

07 

where, 

f x nose 

Ei= x t h(x)dx ,i = 0,1,2, 3, 4 (42) 

J ^tail 

By using this approximation we see that the expansion of the cross flow inte- 
grals give only third order terms, so we have: 

I {2) = / r (2) = 0 (43) 



Now we want to write our matrix equation in state space form, so we have to 
find the inverse of the system matrix: 

Oil £12 Q13 f| 

D D D u 

£21 q 22 023 f| 

D D D u 

031 032 033 p| 5 

D D D u 
0 0 0 1 




where, 



Q-ll — (^22 (I xx Np) ( Ixz ^r)( -^X2 ^Yp) 

12 = (Y^* 'mz g) N p) {Jxz Nr)( rnz G Yp) 

a 13 = {piXQ Yf)( Ixz Np) [I zz Nr)( TYIZq Yp) 

ci2i = (ttiX(7 Ny)(^I X x Np) + ( ttlzq Ky)(^ I xz Np) 
CL22 = (rn — Yy)(Ixx - Kp) - (- mz G ~ Ky)(-mz G ~ Yp) 

&23 — {rn Yy){^ Ixz Np) -f - (rnx G Ny){^ ttizq Ip) 
asi = {mx G ~ Ny)(I xz - Kf) - (- mz G - Ky)(I zz - Nr) 
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o- 32 



= -(m - Yi)(I xz - Kr ) + (- mzc - Ky)(mx G - Yf) 
033 = {m -Yi,){I zz - N+) - (mx G - Ni,)(mx G -Y-r) 

and 

£> = (m - y*)(/« - Nf)(Ixx - K p ) 

- (m - Yi)(I xz - K,){-I xz - Np) 

- (mx G - Nii){mx G - Y+){I XX - Kp) 

+ (mx G - Np)(I xz - Kr){-mz G - Yp) 

+ ( -mz G - Kp)(mx G - Yp)(-I xz - Np) 

- (~mz G - Kp)(I zz - Nf)(-mz G - Yp) 

If we multiply equation ( 31 ) by (A 7 ) -1 from the left side we get: 



(A'y'A'i = (. A')- l B'x + (A')- l g(x) 



or 



x — Fx + G(x) 



where, 



F = (A')~ X B' , 

G(x) = {A')~ Y 9 {x) 

Then matrix F is a 4 x 4 matrix defined as follows: 

“ ^ii ^12 Fv\ Fi 4 " 

1*21 1*2 2 1^3 1*24 

F = £x £2. £2. E± ’ 

D D D D 

0 0 10 



(44) 



(45) 



(46) 

(47) 
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where, 



F n - 

F 12 = 

F 13 = 

Fu = 

F 21 = 

F22 = 

F 23 = 

F 24 = 

-F 3 i — 
■P32 = 
F 33 = 
-F34 = 

From equations 



where, 



a n (Y v U) + a 12 {N v U ) + a 13 (/^t/) 

aii(F r — m)U + a 12 (N r — mx G )U + a i 3 (K r + mz G )U 

a n{Y p U) + a\ 2 {N p U) 4 - ffli3 {KpU) 

an (IF - -6) + ai2(xG^F - xbB ) + <ii 3 (-.zgW + zbB ) 

a 2 i(y„C/) + a 22 (N v U ) + a 23 (K v U) 

«2i(F r — m)U + a 22 ( N r - mx G )U + a 23 ( K T 4 - mz G )U 
a 2 i ( Y p U ) + a 22 ( N p U ) + a 23 { K p U ) 
a 2 i(VF -5)4- a 22 (xclF - x b B ) + a 23 (-XG]F + z B 5) 
031 (nt/) + a 32 ( N v U ) + a 33 ( K v U ) 

031 (Yr ~ m)U + a 32 (N r - mx G )U 4 - a 33 (K r + mz G )U 
031 (YpU) + a 32 (N p U ) + a 33 (K p U ) 

03i(lF - B) + a 32 {x G W - xbB) + a 33 (-z G W + zbB ) 



( 45 ), ( 47 ) we see that: 



G{x) = ( A ') 1 g(x) = 



G i(x) 
G 2 (x) 
G 3 ( x ) 
G 4 { x ) 



5 



Gi(x) = ^9i(x) + ^g 2 (x) + ?fg 3 (x) 

G2(x) = ^-gi(x) + ^-g2(x) + -j^-gz{.x) 

G 3 ( x ) = ^ gi ( x ) + ? fg 2 ( x ) + ^ g 3 ( x ) 
G 4 (x) = 0 
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Equation (45) can also be written in the following way: 



i = Fx + G { 2 \x) + G { 3 \x) + k (48) 



where, 



and 



G (i >(x) = 



<4 2> w 

G?\x) 

<??>(*) 

G^\x) 




G?(x) 

Gf\x) 

G?\x) 

G?\x) 



) 



) 



k i 




k 4 



Each element of the above non linear terms vectors can be written as follows: 



and 



ofM 

G?\x) 

G?\x) 

af(x) 



ail (2) / v a 12 (2)/ x , ai3 (2)/ x 
~q9 l ( x ) + -Jj-92 (z) + -Jj-93 ( x ) 

°21 (2)/ \ . a 22 (2)/ \ , 023 (2)/ x 

-p-^1 0) + ^52 W + -p93 (z) 

031 (2)/ \ . O 32 (2)/ \ . 033 (2)/ x 

-p"01 (*) + -^-^2 W + ~~^"9z (*) 



^(3)/ \ On (3)/ \ , O 12 (3)/ x , O 13 (3)/ x 

G l \z) = —g x \x) + —g\ '( x ) + —S 3 (*) 

^-,(3) / x 021 (3) / \ , 022 (3)/ x , 0 2 3 (3)/ x 

G 2 (z) = -^9\’(z) + —g\ 1 (x) + —g\ , {x) 

M 3 )f \ 031 (3 ) , x . O32 (3) O33 (3) 

G Y\z) = —g\\x) + —g\ 1 (x) + —g\’(x) 

Gf\x) = 0 
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Finally, the constant terms are: 



h 

k 2 

h 

&4 



an ai2 fli3 

~D Cl + T° 2 + ~D C3 

a 2l . a 22 , ^23 

~5°' + T c 2 + ~d C3 

a 31 , ®32 , 033 

T Cl + ~d 02 + T ca 



C. COORDINATE TRANSFORMATIONS 

From equations (45) and (48) it is obvious that the stability of our system 
depends on the eigenvalues of matrix F. Since we want to investigate the 
behavior of our system oround the Hopf Bifurcation point, it is useful to bring 
our system in its normal coordinate form. This can be done by applying a 
transformation of the coordinate system, using as transformation matrix T, 
the modal matrix of eigenvectors of F, evaluated at a critical point. This 
critical point will be a pair of zq and xq values, that belong to the critical line 
of loss of stability. The applied transformation will be as follows: 

x = Tz (49) 

or, 

^ = T~ l x (50) 

Then equation (48) can be written as follows: 

Ti = FTz + G ( - 2) {Tz) + G {3 \Tz) + k (51) 
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or, 



i = T~ l FTz + T~'G {2) (Tz ) + T~ l G {z \Tz) + T~ l k 



(52) 



In this new coordinate system, the system’s matrix is T l FT , and at the Hopf 
Bifurcation point can be written as follows: 



T -1 FT = 



0 -UQ 0 0 

0 0 0 
0 0 Pi 0 

0 0 0 P 2 



Where uq is the imaginary part of the critical pair of eigenvalues, and the 
remaining eigenvalues Pi , P 2 are negative. For values close to the Hopf Bi- 
furcation line, we can write the system matrix as follows: 



T -1 FT = 



a! t 



— (wo "h OJ* e) 0 0 

(u>o + w'e) a'e 0 0 

0 0 (Pi + Pj'e) 0 

0 0 0 (P 2 + P^e) 



where: e = Difference from the critical point (zq — zq c ). 



a' = Derivative of the real part of the critical eigenvalue with respect to e. 

u' = Derivative of the imaginary part of the critical eigenvaluewith respect 
to e. 



P[ = Derivative of Pi with respect to e. 
P' 2 = Derivative of P 2 with respect to e. 



D. CENTER MANIFOLD EXPANSIONS 

By using the coordinate transformation described in the previous chapter, 
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and from equation (50), we see that we have a new variables vector, which is: 

’ zi ' 

z= 22 

z 3 

. 2 4 . 

The first two coordinates z\ , 22 , are the critical coordinates and correspond 
to a pair of complex conjugate eigenvalues. The remaining two coordinates 
23 , z 4 , are the stable coordinates and they always correspond to eigenvalues 
that are negative. The center manifold theory predicts that the relationship 
between the critical coordinates z\, 22 , and the stable coordinates 23 , 24 , is at 
least of quadratic order. After this assumption, the two stable coordinates can 
be written as follows: 



Z3 = h n z\ + h 12 2 i 2 2 + h 22 zl 


(53) 


2 4 = Sn2 2 + S 12 2i2 2 + S22z\ 


(54) 



The coefficients hij, Sij, of equations (53), (54), need to be determined. By 
differentiating equations (53), (54), we get the following: 

23 = 2h\\Z\Zi + hi 2 (ziZ 2 + 2122 ) + 2 J 122 Z 2 Z 2 (55) 

24 — 2snZ\Zi + S\ 2 {z,\Z 2 + Z 1 Z 2 ) + 2S22 Z 2 Z 2 (56) 

From equation (52) and if we ignore the higher order terms, we take: 



Z\ = <^0 z 2 


(57) 


2 2 = l^OZl 


(58) 



If we substitute equations (57), (58), into equations (55), (56), we get the 
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following: 

1 3 = hx^jjQzl + 2(/i 2 2 - h u )ujQZiZ2 - huUJQzj (59) 

1 4 = S 12 UJ 0 Z 1 + 2(s 2 2 - su)uQZiz 2 - s\ 2 ujoz% (60) 

The third and the fourth equations of matrix equation (52) can be written in 
the following way: 

is = Pxzt + lT-'GWiTz^+T-'kw (61) 

Z4 = P2Z 4 + [T- l G^(Tz)] {4A) +T- 1 k {4A) (62) 

In equations (61), (62), we kept terms up to second order. 

The transformation matrix T and it’s inverse T~ l are 4x4 matrices, and 
their elements can be denoted by: 

T = [my], T _1 = [ny], i, j = 1,2, 3, 4 (63) 

From equation (52) we have: 

d\ 

d 3 

d A 

where, 

di = n u G { ?\Tz) + nnG%\Tz) + n u Gf\Tz) 
d 2 = n 21 G { x\Tz) + n 22 G { 2 \Tz) + n 23 Gf\Tz) 
dz = n zl G { i\Tz) + n z2 G ( i\Tz) + n^Gf\Tz) 
d A = n A xG [ x\Tz) + n A2 G { 2 \Tz) + n A zG ( i\Tz) 



(64) 

(65) 

( 66 ) 
(67) 



T~ l G {2) {Tz) = 
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Also from coordinate transformation, the relationship between the old and 
new coordinates is as follows: 



v — mnzi + mi 2 Z 2 + m 13 z 3 + 7ni 4 z 4 (68) 

r = m 2 i2i + 7n 22 z 2 + m 23 z 3 + tu 24 z 4 (69) 

p = mzizi + m 32 z 2 + m 33 z 3 4- m 34 z 4 (70) 

<f> = m 4i zi + m 42 z 2 + m 43 z 3 + m 44 z 4 (71) 



Finally if we substitute equations (68), (69), (70), (71), and the expressions for 
G i, G 2 , G 3 , G 4 into equations (64), (65), (66), (67), we get the final expressions 
for the coefficients 

d\ — rcn(Zi5Zj + Zi6-2i-22 + h7 z \) 

+ ni2{h$ z \ + ^26 z i z 2 + ^27-2 2 ) 

+ n \z{hb z \ + he z i z 2 + ^37-2 2 ) (72) 

d-2. = n 2l(h5 z i + h6 z \ z 2 + h7 z 2) 

+ 7722 ^ 25-21 + hz z l z 2 + ^ 27 - 2 2 ) 

+ 7 l2z{hb z \ + he> z \ z 2 + ^37-2 2 ) (73) 

dz = 7z 3 i(/i 5 zi + Zi6-zi-Z2 + ^172 2 ) 

+ 7Z 32 (Z 2 5Z2 + hz z \ z 2 + ^27-2 2 ) 

+ 77 33 (Z 3 5Zj + h6 z l z 2 + h7 z %) (74) 

d 4 = 77 4 i(Zi5Zi + l l6-2lZ 2 + l\7 z 2) 

+ 77 42 (/ 25 Zi + l2Z z \ z 2 + h7 z 2 ) 

+ ti 43 (/ 3 5Z2 + he z i z 2 + h7 z 2.) (75) 
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Coefficients kj,i — 1,2,3 j — 5,6,7 in equations (72), (73), (74), (75), are as 
follows: 

^1,5 = -^-yc(ra 3 l + m 2l) 

+ ~J^{Ixy m 3i + 7 y2 m 3im2i + J/G m 21 m ll) 

- ^(/i y m 3 im2i + /y^m^ + my G miim 3 i 

+ ~ VbBttiI^ (76) 

^ 1,6 = ~ 2 /G( 2 m 3 im 3 2 + 2m 2 i + m 2 2 ) 

+ ~“[2/i y Tn 3 im 32 + I yz {mz im 22 + m 32 m 2 1 ) 

+ yG(^2i^i2 + m 22 mn)] 

0,-^3 

- — [/xy(^3l"2-22 + m 32 m 2l) + 2 Jy 2 777 2 l77722 

+ my G (mum 32 + mi 2 m 3 i) + 2(y G W - yB£)m 4 im 42 ] (77) 

, a n t 2 . 2 \ 

*1,7 = -^-2/<?(™ 32 + ^ 22 ) 

+ ^(/xy7n3 2 + /y 2 m 32 m 22 + j/ G m 22 7n 12 ) 

- ^[7xym 32 m 22 + 7yx 772-22 + 7777/(3771 12777 32 + (y^VK - y B . 8 ) 77742 ] ( 78 ) 

. a 21 / 2 , 2 \ 

*2,5 = ^-2/g(777 31 + 777 21 ) 

+ ^(/xy7773! + /y 2 777 3 i777 2 l + 1/0^21^111 ) 

- — (7xy777 3 l7T7 2 l + /y 2 777 2 i + 777y G 777 n 777 3 i 

+ yoH 7 777^ - ysB 777^) ( 79 ) 

*2,6 = ~^Vg( 277731777 3 2 + 277721 + 7772 2 ) 

+ --^-[2/xy7773i777 3 2 + -*yz (77731 777 22 + 77732777 2 l) 

+ yG(7n 2 l777i2 + 777 2 2777 ii)] 
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- -Jj[lxy{ m Zl m 22 + 77l 3 277i 2 i) + 2Iyzm2imv2 

+ myG(mnm 3 2 + mi 2 m 3 i) + 2(y G W - 7 / 5 . 8 ) 7714177142 ] (80) 

*2,7 — ~j^yG\ m 32 + m 22) 

+ ~JJ~{Ixy m 32 + IyzT n 32 m 2‘2 + 7/G771 2277112) 

- -^[4^32^22 + IyzTn.22 + rnycm^m^ + (vgW - yBB)m\ 2 ] (81) 

*3,5 = jfVGimh + m|j) 

+ -^-(Ixyml 1 + /y 2 m 3 im2i + 2/G^21^1l) 

- -^■(Ix y msim2i + Iy Z m 21 + mt/G"iii"i3i 

+ yciy - yB^m^) (82) 

*3,6 = ^-yG(2m 3 im 32 + 2m 2 i + m 22 ) 

+ ^[27x^7713177132 + 7y 2 (77131771 22 + 771 32 77l2l) 

+ 2/g(" 121^12 + 77i2277ill)] 

- ^[^(7713177122 + 771 3 277l2l) + 27 yz 77l 2 l77l22 

+ my G (m n m 3 2 + 77 ii 2 77 i 3 i) + 2(y G W - 7 / 58 ) 7714177142 ] (83) 

*3 , 7 = ^7/G(77l32 + 77l2 2 ) 

+ ^(7 ly 77l3 2 + /j, 2 771 3 277l22 +yG77l22771l2) 

- -“[-*12/7713277122 + -fj/zTTl^ + 7717/G771i277i32 + (l /G^ ~ 7 / 5 . 8 ) 77142 ] (84) 
Then equations (61), (62), can be written as follows: 



^3 — -P1Z3 + 7*3 + e 3 



— -82-24 +7*4 + 64 
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( 85 ) 

(86) 



where e 3, come from the constant terms as: 



and 



T~ l k = 



e\ 

e2 

S3 

e4 



e\ = k\7l\\ + ^2^12 + ^37113 

e 2 = ^l n 21 + ^2^22 + k 3 Tl2Z 

e 3 = fcin 3 i + k^nz 2 + k 3 n 3 z 

e 4 = kin 4 i + k 2 n 4 2 + k 3 n 43 

Using equations (53), (54), we can write equations (85), (86), as follows: 

£3 = P\{h\\z\ + h\ 2 z\Z2 + h 2 2Z 2 ) + d 3 e 3 ( 87 ) 

Z 4 = ^(Sll^l 4" S\2Z\Z 2 + S 22 Z 2 ) + 6/4 + 64 (88) 

And using equations (72), (73), (74), (75), they become: 

£3 = (Plhn + U31Z15 + 7132^25 + n 33^35) 2 l 

+ (Plh\2 + 7*31^16 + n 32^26 + ^33^36)^1^2 

+ (-^ 1^22 + n 3 Ul 7 + n 32^27 + 7133/37)22 + 63 ( 89 ) 

£4 = (P2S1I + ^ 31^15 + n Z 2^25 + 72.33/35)2^ 

+ (^2512 + 77.3 i / i 6 + 77.32/26 + 7733^36)^122 

+ {P2S22 + 7731/17 + 7732/27 + 7733/37)2! + 64 ( 90 ) 

Comparing the coefficients of equations (89), (90), with the coefficients of 

equations (59), (60), we get: 

— P\hu + u>()h 12 = 7731/15 + 7732/25 + 7733/35 
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— 2 u>ohn — P \ h \2 + 2 uJoh 22 = ^ 31^16 + 7732/26 + W33/36 
— UQ)h 12 — P 1/122 = n Z 1^17 + W32/27 + 7733/37 

Solution of the above 3x3 linear system of equations, gives us coefficients /i n , 

h \ 2 , h 2 2 - 

Also in the same way: 

—P2S11 + W0S12 = 72.41/15 + 7742/25 + 7743/35 
— 2 u;oSii — P2S12 + 2 u!qS 22 = 7741/16 + 7742/26 + 7743/36 
— o;o 5 i 2 — P2S22 = 7741/17 + 7742/27 + 7743/37 

Solution of the above 3x3 linear system of equations, gives us coefficients sn, 
Sl2) s 22 - 



E. AVERAGING 



In this part of our analysis we are going to take into account the third order 
terms of equation (52): 



T - 1 G ( 3 ) (Tz) 



fl 

h 
h 
f* . 



Where, 



fi — nnG^(Tz) + ni2G2(3)(Tz) + n\zG^\T z) (91) 

h = 772i<?i 3) (Tz) + n 2 2G 2 (3)(Tz) + n 23 G$\Tz) (92) 

h = 773 lGS 3) (Tz)+ 7732 G 2 ( 3 )(rz)+ 7733 Gi 3 ) (Tz) (93) 
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h = n 4l Gf\Tz)+n 42 G 2 (3)(Tz)+n 43 G^(T:z) 



(3), 



(94) 



(95) 



(96) 



If we substitute equations (68), (69), (70), (71) and the expressions for G i, 
G 2, G3 into equations (91) throuh (94), we get: 

/1 = ra ii(^ii- 2 i + h 2 z\z 2 + l\z z \ z 2 + ^14^2) 

+ ^12(^21^1 + h 2 z \ z 2 + hz z \ z 2 + /24 ^2 ) 

+ n iz(hi z i + h 2 z i z 2 + hz z i z 2 + ^34 z 2 ) 

/2 = ra 2l(^ll 2 i + ^12'2 i22 + hz z \ z 2 + ll4 z 2 ) 

+ 7Z22(^21 2 i + h2 z i z 2 + hz z \ z 2 + h4 z 2 ) 

+ ^23(^31^1 + h 2 z \ z 2 + hz z l z 2 + 134-22) 
fz = ra 3l(^ll-2l + h 2 z \ z 2 + hz z \ z 2 + h4 z 2 ) 

+ n Z 2 ^ 21 z 4 + l 22 z \ z 2 + hz z l z 2 + ^24 z 2 ) 

+ ^33(^3l2i + lz 2 z \ z 2 + hz z \ z 2 + h4 z 2 ) 
f 4 — n 4 \{l\\z\ + l\ 2 z \ z 2 + hz z l z 2 + h4 z 2.) 

+ n 42 {l 2 \z\ + l22 z i z 2 + hz z l z 2 + h4 z 2 ) 

+ ^43(^31 2 i + h 2 z i z 2 + hz z \ z 2 + h4 z l) 

From equation (52) and from the system matrix, the derivatives of the first 
two modified coordinates, i\ and Z 2 can be written as follows: 



(97) 



(98) 



z\ = {a'e)z\ — (cao + u>'e)z 2 + FF\(z\, z 2 ) 
z 2 = ( w o + u)' e)z\ + a! ez 2 + FF 2 (z\, z 2 ) 



(99) 

( 100 ) 



where, 



F F i(z\, z 2 ) = di + h + ei 



( 101 ) 
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FF 2 {zi, 22 ) = d 2 + h + e 2 



( 102 ) 



If we combine equations (101) and (102) with equations (72) through (75), 
and (95) through (98), we get: 

F Fi(zi, z 2 ) = r u zl + r 12 zlz 2 + r iz zizl + ruzl 

+ Pn z i + VV2 Z \Z 2 + P 13 Z 2 + e\ (103) 

FF 2 (z 1 , z 2 ) = r 21 z\ + r 22 z\z 2 + r<az\z\ + r 2 A z\ 

+ P21*1 + P22Z\Z 2 + p 2 zzl + ei (104) 

where coefficients and pij are: 

= ^11^11 + 7212/21 + 7213/31 
ri 2 = 72.11/12 + 7212/22 + 7213/32 
r 13 = 72ll/l3 + 7212/23 + 7213/33 
r 14 = 72 11 /14 + 7212/24 + 7213/34 
^21 = 72 2l/ll + 7222^21 + 7223/31 
7" 22 = 722i/i2 + 7222^22 + 7223^32 
7*23 = 7221^13 + 7222^23 + 2223/33 
r 24 = 7221^14 + 7222^24 + 2223/34 

or generally, 

7* ij = nnhj + ni 2 l 2 j + n^l^j i = 1,2 j = 1,2, 3,4 (105) 

also, 

Pn = 72 n/i 5 + 2212/25 + 2213/35 



39 



Pl2 ~ nilll6 + ^12^26 + ^13^36 
Pl3 = n nh7 + 7712^27 + 7113^37 
P21 = 72-21^15 + 7122^25 + 7723^35 
P22 ~ 72-21^16 + 72.22^26 + 72-23^36 
P23 = 7121^17 + 7122^27 + 72.23^37 



or generally, 



Ptj = Tiji^u + nuhk + 7ii3/3fc 7 — 1,2 j = 1,2,3 k=j+ 4 



(106) 



The coefficients Uj i = 1, 2, 3 j = 1, 2, 3, 4 are as follows: 

^li = — zr\ z ~ (E 0 m 3 n + 3^1777^77721 + 3E 2 m n m 2 21 + Ezm \ x ) 
L ) 07 

— -— [ + 3B2 m u m 2l + 3Esmum2i + £^7, ) 

L7 07 



1 



+ ^G^ - XBB)Tn z Al ] + ^p-(z G W - z B B)m A1 

0 0 U 

li2 = --^■[-^■(SE 0 m 2 n m l 2 + 3^1(777^77722 + 2 m n mi 2 m 2 i) 

u 07 

1 



(107) 



+ 3 ^ 2(7771277721 + 2m 21777 22777 ii) + 3£3m 21 m 22 ) + -(W — B)3m A1 m4 2 ] 

6 

- ^[^T L (3Eim 2 ll m 12 + 3E 2 ({Tn\ l m 22 + 2m n mi 2 m 21 ) 

V 07 

+ 3£3((m 12 TO2 1 + 2m 2 im 2 2777ii) + 3 E^m^m^) 



1 



Ol3 



+ ~(x G W - x B B)3m4 1 m4 2 ] + — (z G W - z B B)3mi 1 m 42 



6 



6 D 



(108) 



/13 = -^-[^ L (3 J Eo777iimi 2 + 3£; 1 (mi 2 m 2 i + 2 mnmi 2 m 2 2 ) 

JJ 07 

+ 3£2(777 22 m u + 2m 2 im 22 mi 2 ) + 3Ezm 2 \m\ 2 ) 



+ ^(W 7 - B)3m 4 iml 2 } 
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- 3£?iJnnmf 2 + ZE2(rn\ 2 m 2 i + 2m 11 mi 2 m 2 2) 

07 

+ 3i?3(m 22 mii + 2m 2 im 22 mi 2 ) + 3f?4m 2 im| 2 ) 



1 



+ - rcBB)3m 4 im4 2 ] + ^{z G W - z B B)3m 4l m 2 42 

S~1 

li4 = — FT ( E ° m 12 + 3£ i m n m 2 i + 3£ 2 7n 12 m 22 + £ 3 m 22 ) 

07 

+ ^(W-B)m y 

- ^[^£2.(£jmf 2 + 3 F 2 ™?i>H 2 i + 3 Esm 12^22 + £ 4 ™ 2 2 ) 

07 

+ i(x G ^ - x B B)m\ 2 ) + ^(z G H/ - z B B)m\ 2 



021 r C_Dy 



6D 



hi - — — [— — (£ 010 ^ + 3^1771^77121 + 3£ 2 ioiim 21 + £ 3 m 21 ) 

V D7 

+ i(IV-B)m^] 

- -^rl-TT-MBarctii + 3 B 2 mfim 2 i + 3 B 3 m n m 21 + £ 4 ro 21 ) 

07 

+ ^ {xgW - x b B + ^(z G W - z B B)m 3 41 

*22 = -^r[^-{^E 0 m 2 n m 1 2 + 3Ei(m 2 n m 22 + 2 m 11 mi 2 m 2 i) 

07 



(109) 



( 110 ) 



(111) 



+ 3^2(1012^21 + 2m21m22io.11) + 3 £J 3 m 21 m 2 2 ) + 7(1^ — B)3m 2 4l m 4 2\ 

0 

_ \9^(^,E l m\ l mi2 + 3£^ 2 ((iOi 1 m22 + 2m 1 im 12 m 2 i) 

L) 07 

+ 3^ 3 ((mi2m 21 + 2m 2 im22iOii) + 3^4m 21 m 2 2) 



1 



023 



+ t;(xgW - x B B)3m 2 4l m 42 ] + - z B B)3m z 4l m 4 2 

0 oL> 

*23 = -^■[^^-(3Eom u ml 2 + 3E 1 (m1 2 m2i + 2m n mi2m22) 

JJ 07 

+ 3 ^ 2(10221011 + 2m2im 2 2ioi2) + 3Ezm2\m 22 ) 



( 112 ) 



+ -{W - B)3m 4 im 2 42 ] 

6 

— -^-[-^-{3Eim n m 2 2 + 3E2(rn\ 2 rri2i + 2m nmi 2 m22) 

D 07 
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+ 3Ez(ml 2 mn + Zrrizimyzmu) + 



+ ^(xgW - x B B)3m 4 iml 2 } + ^|( z G W - z B B)3m 44 m\ 2 
o o V 

1 24 = -^[-^(Eomi2 + 3Eim\ l m2\ + 3 ^ 2 ^ 12^22 + Bzm\ 2 ) 



(113) 



D 67 



1 



+ |( W-B)m\ 2 ) 

- -^\-^-{Eim\ 2 + 3^2^1477721 + 3^37771277722 + ^ 4 ^ 22 ) 

V 07 

+ ifioH 7 - x B B)m 3 a2 ] + ^L(z a W - z B B)m \ 2 



(114) 



/ 3 I = — [ y ■)“ 3 £?i777 14 77721 + 3 £^2 m 'll' m '21 ■£'3 r 77 ’2l) 



D 67 



1 



+ |(^-5)mJi] 

a32 ’[“7^(-El m ll + 3^2^11^721 + 3^37771177721 + £^4"72l) 



£> 67 

1 



+ E( XgW _ o; Bj B)777^i] + ^r(z G W - z b B ) 



m 
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(115) 



6'”~ " ' '« 1J ' 6£> 

^32 = ( 3 E 0 m 2 n mi 2 + 3 £i (777^77722 + 27 nur 77 i 27 n 2 i) 

D 07 

+ 3£^2(”7l2?7721 + 2?77 21777 22777n) + 3£37772i77722) + g(^ — B)3m\ l m 4 -^ 

- ^[^-(3£:i777ii777i2 + 3£ 2 ((777 ii 777 2 2 + 2777n 777 i2777 2 l) 

V 07 



+ 3^3 (( 777 1277721 + 2777217772277711 ) + 3£ 4 7772i777 2 2) 

+ \{. x gW - x B B)3m 2 4l m 4 ^ + ^(zqW - z B B)3m 2 4l m 4 2 

O vV 

kt - -^■[-^■(3^07771177712 + 3£ , i(777 i 2 77721 + 2777117771277722) 

+ 3 £/2(777 22 777 11 + 2777217772277712) + 3Ez'm2\'m\ 2 ) 

+ ^(M 7 - B)3777 4 i777 42 ] 

_ ^]^^^Eirn U m\ 2 + 3£;2(777i 2 77721 + 2777n 777 i 2 77722) 

V 07 

+ 3 £' 3(77722 777 11 + 2777217772277712) + 3£ 4 7772l7772 2 ) 



(116) 
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(117) 



hA 



+ \{. x gW ~ x B B)3m 41 ml 2 ] + ^:{.z G W - z B B)3m 4 im\ 2 

o oV 

_ _^.[^i(£: om J 2 + 3£? 1 mi 1 m 2 i + 3£ 2 mi2ra22 + E 3 ml 2 ) 

V 07 

+ \{W-B)m\z) 

- ■7rlir 2 '(£i" 1 i2 + 3£2m;,m2i + 3E 3 m 32 ml 2 + E A m! n ) 

D 07 



+ ^{x G W - x B B)m \ 2 ] + ■jr^( z GW ~ z B B)m \ 2 
The next step is to introduce polar coordinates in the form: 



( 118 ) 



z\ = R cos# (119) 

Z2 = -Rsin0 (120) 

We use polar coordinates, because it is easier in this way to investigate the 
existence of limit cycles. 

Substituting equations (119), (120), into equations (99), (100), we get: 

R cos 8 — i?(sin 9)9 = (wo + a/e)i?cos0 + t/e-Rsin# + 

P 4 (9)R 3 + Q 1 (9)R 2 + e 4 (121) 

R sin 9 + R(cos9)9 = (cuo + u'e)R cos 8 + a eR sin 9 + 

P 2 ( 0 )R 3 + Q 2 ( 0 )R‘ 2 + e 2 (122) 

where, 

Pi(9) = rn cos 3 9 + t~i2 cos 2 9 sin 9 + r 2 3 cos 9 sin 2 9 + r 44 sin 3 8 (123) 

Q\{9) — P 11 cos 2 0 +pi2cos0sin# + pi 3 sin 2 6> (124) 

P 2 (0) = r -21 cos 3 9 + r 22 cos 2 9 sin 9 + r^z cos 8 sin 2 9 + r 24 sin 3 9 (125) 

QziQ) = P 21 cos 2 9 + P 22 cos 9 sin 9 + p 2 3 sin 2 8 (126) 
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If we multiply equation (121) by cos# and equation (122) by sin#, and add 
the two resulting equations, we get: 

R = a'eR 4- P{6)R 3 + Q(9)R 2 + (ei cos # + e 2 sin #) (127) 

where, 

P{9) = Pi(#) cos# + P 2 (#)sin# (128) 

Q(0) = Qi(#)cos# + Q 2 {9) sin 9 (129) 

Equation (127) contains one variable that varies slowly in time ( R ) and a fast 
variable (#). 

If we average this equation over one complete cycle in #, from 0 to 2-7T, 
equation (127) becomes: 

R = a'eR + KR 3 + LR 2 4- M (130) 

where, 

* - 



— 0 (3r*n 4- + r 22 + 3r24 


(131) 


1 


r27T 




L = — 

27 r j 


/ Q(#)d# 

'0 


(132) 


1 






m - s. 


/ ei cos # d# 




1 


r27r 




27T J 


/ e 2 sin # dO 
'0 





= |^[sin#]^ - ^[cos#]^ 



= 0 (133) 
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Finally equation (135) becomes: 



R = aeR + K R 3 



(134) 



F. LIMIT CYCLE ANALYSIS 

At steady state R = 0,and equation (134) becomes: 

0 = R(ae + KR 2 ) (135) 

Equation (135) has two solutions. The first solution is R = 0. This is the 
trivial solution and it does not give us much information. The second solution 
is: 

R = ^ (136) 

This solution gives us a limit cycle of constant amplitude R in the z\, Z 2 
cartesian coordinate system. This limit cycle exists if the quantity inside the 
square root is possitive, or 

> 0 (137) 

Condition (137) is necessary for the amplitude of the limit cycle, R, to be a 
real number. 

In our case a' is always negative, because for constant as we decrease 
e (Figure 1), the real part of the critical pair of eigenvalues increases, due to 
further loss of stability. In other words we can say that: 

a < 0 (138) 
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From conditions (137), (138) we see that the existence of the limit cycles 
depends on the value of parameter K. We can see that: 

• If K <0, periodic solutions exist for e < 0 or zq — zq c < 0 or zq < zg c , 
and they are stable. 

• If K > 0, periodic solutions exist for e > 0 or zq — zg c > 0 or zq > zq c , 
and they are unstable. 

The characteristic root of equation (134) in the vicinity of (136) is: 

/3 = -2 at (139) 

The sign of this characteristic root, assigns the stability of the periodic solu- 
tions. 



G. RESULTS AND DISCUSSION 

Typical results in terms of the nonlinear stability coefficient K are pre- 
sented in Figures 3 through 6. The stability coefficient K is shown in its 
normalized form as K ■ 7 (Papadimitriou, 1994). Figure 3 shows K versus 
the LCG/LCB separation distance xq (in ft) for a given vehicle speed and for 
different values of the drag coefficient. It can be seen that K is everywhere 
negative, which means that all bifurcations to periodic solutions are super- 
critical. Higher values of the drag coefficient result in stronger supercritical 
bifurcations which means that the corresponding limit cycle amplitudes will 
be smaller. 
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K'GAMMA K’GAMMA 




Figure 3: K ■ 7 versus xq for U = 5 ft/sec and different values of Cjj. / 




Figure 4: K • 7 versus xq for = 0.5 and different values of U (ft/sec) 
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Figure 5: Simulation results (</>,£) for C]j y = 0.5, U = 5 ft/sec, and xq — 1 ft 




Figure 6: Limit cycle amplitudes (in <t>) versus zq for U = 5 ft/sec and xq = 1 ft 
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Figure 4 shows K versus xg for a given drag coefficient and for different 
forward speeds. It can be seen that the bifurcations are supercritical with the 
possible exception of large speeds and small xg. This means that it is possible 
for a properly trimmed vehicle at relatively high speeds to experience an oscil- 
latory behavior even before stability is lost. This demonstrates a destabilizing 
effect which could not have been predicted by linear techniques. 

Figures 5 and 6 shows numerical simulation results which confirm the theo- 
retical predictions. Both figures correspond to a supercritical bifurcation case 
and they show a continuous increase in limit cycle amplitudes as the bifurca- 
tion point is crossed. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



This thesis presented a comprehensive nonlinear study of straight line sta- 
bility of motion of submersibles in coupled sway/yaw/roll motions under open 
loop conditions. Primary loss of stability was shown to occur in the form of 
Hopf bifurcations to periodic solutions. This loss of stability is characteristic 
of the coupling of roll into sway and yaw and cannot be predicted by consid- 
ering the uncoupled motions. The critical point where instability occurs was 
computed in terms of vehicle metacentric height, longitudinal separation of 
the centers of buoyancy and gravity, and the forward speed. Analysis of the 
periodic solutions that resulted from the Hopf bifurcations was accomplished 
through Taylor expansions, up to third order, of the equations of motion. A 
consistent approximation, utilizing the generalized gradient, was used to study 
the non-analytic quadratic cross flow integral drag terms. The results indi- 
cated that loss of stability occurs always in the form of supercritical Hopf 
bifurcations with stable limit cycles. It was shown that this is mainly due to 
the stabilizing effect of the drag forces at high angles of attack. Sub critical 
bifurcations, however, with considerably higher limit cycle amplitudes may 
develop for sufficiently high forward speeds and small LCG/LCB separations. 
Simulation studies of these sub critical bifurcations along with the effects of 
vertical plane coupling constitute recommendations for further research. 
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APPENDIX 



The following is a list and description of the computer programs used in this 
thesis. The programs are written in FORTRAN or MATLAB. Complete print- 
outs of the programs follow after the list. 

• STABILITY. M 

MATLAB program for performing linear stability analysis. 

• SIM.M 

MATLAB program and functions for numerical simulation. 

• HOPF.FOR 

FORTRAN program for calculating the nonlinear stability analysis coef- 
ficient K. It requires data from STABILITY.M and standard numerical 
linear algebra subroutines. 
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7. STABILITY. M 

7. 

7. LOSS OF STABILITY 

% ***************************************** 
a=l 

W=12000; 

IXX=1760 ; IYY=9450 ; 

IZZ=10700 ; IXZ=0 ; IXY=0; 

IYZ=0 ; L=17.425; RH0=1.94; 

G=32 . 2 ; U=6 . 5 ; M=U/G; B=W; 

ND1=0.5*RH0*L'2; 

7. DEFINE HYDRODYNAMIC COEFFICIENTS 

YPDOT=l . 270e-04*NDl*L“2 ; 

YVD0T=-5 . 550e-02*NDl*L ; 

YRDOT=l . 240e-03*NDl*L~2 ; 

YP=3 . 055e-03*NDl*L ; 

YV=-9 . 310e-02*NDl ; 

YR=-5.940e-02*NDl*L; 

NPD0T=-3 . 370e-05*NDl*L"3; 

NVDOT=l . 240e-03*NDl*L~2; 

NRDOT=-3 . 400e-03*NDl*L~3; 

NP=-8 . 405e-04*NDl*L‘2 ; 

NV=-1 . 484e-02*NDl*L ; 

NR=-1 . 640e-02*NDl*L"2 ; 

KPDOT=-l .01e-03*NDl*L"3; 

KVDOT=l . 27e-04*NDl*L‘2 ; 

KRDOT=-3 . 37e-05*NDl*L"3; 

KP=-1 . 10e-02*NDl*L‘2 ; 

KV=3 . 055e-03*NDl*L ; 

KR=-8 .41e-04*NDl*L‘2 ; 

flag=0; 

for XG=0 : 0 . 01 : 2 , 
f lag=f lag+1 ; 
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xg(f lag)=XG; 

a=IXX-KPDOT; b=KP*U; e=KV*U; 
f =KRDOT ; i=YP*U; j =M-YVDOT ; k=YV*U; 

1=XG*M-YRD0T ; m=U*(YR-M); o=NPDOT; 
p=NP*U; q=-XG*W ; r=XG*M-NVDOT; 
w=U* (NR-XG*M) ; x=NV*U; u=IZZ-NRDOT; 

al=-u*M~2 ; 

a2=-u*M*YPD0T-u*M*KVD0T+M*o*l+f *r*M; 

a3=a*j*u-a*l*r-u*KVD0T*YPD0T+KVD0T*o*l+f *r*YPDOT-f *o*j ; 
bl=(w* (M~2) )+(r*U* (M~2) ) ; 

b2=-M*e*u-M*u*i+w*M*YPD0T+w*KVD0T*M-o*m*M+p*l*M-f *x*M+. . . 
r*M*U*YPDOT-o* j *M*U+r*KR*U*M ; 

b3=-e*u*YPD0T+e*o*l-u*i*KVD0T+w*KVD0T*YPD0T-KVD0T*o*m+KVD0T*p*l 
a*j*w-a*k*u+a*l*x+a*r*m-b*j*u+b*l*r+f*r*i-f*x*YPDOT+. . . 
o*k*f-f *p* j+r*KR*U*YPDOT ; 

cl=-x* (M~2) *U; 

c2=r*i*M*U-x*M*U*YPD0T-x*KR*U*M+o*k*M*U-p* j *M*U+ j *u*W-l*r*W- . . . 
q*l*M+w*i*M-m*p*M+e*w*M ; 

c3=a*k*w-a*m*x+b* j *w+b*k*u-b*l*x-b*r*m+r*i*KR*U-x*KR*U*YPDOT+. . 
o*k*KR*U-p* j*KR*U-q*l*KVDOT+w*i*KVDOT-m*p*KVDOT-e*u*i+e*w*YPDOT 
e*o*m+e*p*l+f*q* j-f *x*i+f*p*k; 

d2=q* j *M*U-x*i*M*U+p*k*M*U+q*m*M-j *w*W-k*u*W+l*x*W+r*m*W ; 
d3=q*j*KR*U-x*i*KR*U+p*k*KR*U-f *q*k+q*m*KVDOT-e*q*l+e*w*i- . . . 
e*m*p-b*k*w+b*m*x ; 



e2=k*w*W-m*x*W-q*k*M*U ; 
e3=e*q*m-q*k*KR*U ; 

f 5=(-e2* (bl~2) ) +bl*cl*d2 ; 

f4=( (b2*cl+bl*c2) *d2)+bl*cl*d3-al* (d2~2) -e3* (bl~2) -2*e2*bl*b2 ; 
f 3=-a2* (d2~2)-2*d3*d2*al-2*e3*bl*b2-e2* ( (b2~2)+2*bl*b3)+ . . . 
d2*(b3*cl+b2*c2+bl*c3)+d3*(b2*cl+bl*c2) ; 
f 2=-e3* ( (b2~2)+2*bl*b3) -2*e2*b2*b3+d2* (b3*c2+b2*c3)+ . . . 
d3* (b3*cl+b2*c2+bl*c3)-a3* (d2~2)-2*d3*d2*a2-al* (d3~2) ; 
f I=b3*c3*d2+d3* (b3*c2+b2*c3)-2*e3*b2*b3-e2*(b3~2)-. . . 
2*d3*d2*a3-a2* (d3~2) ; 
f 0=b3*c3*d3-a3* (d3~2) -e3* (b3~2) ; 
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coef = [f 5 f 4 f 3 f 2 fl fO] ; 

ZG=roots (coef) ; 

tot(f lag)=ZG(5, 1) ; 
end 

plot (xg, tot) .grid; 

title (’ZG versus XG plot in the point of loss of stability') 
xlabelC’XG in ft’); 
ylabelC ’ ZG in ft’); 
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7. NON LINEAR SIMULATION PROGRAM 



t0=0 ; 

tf inal=500; 
q=( 1/180) *pi; 
yO= [0 0 0 q] ; 

[t , y] =ode45 ( ’ vdpo2 ’ , tO , tf inal , yO) ; 

figure(l) ,plot (t ( : ,1) , 57 . 29578*y( : ,4)) ,grid 
title (’plot of fi with time’); 
ylabel(’fi’) ,xlabel(’time’) 



57 



% NON LINEAR SIMULATION PROGRAM 



function yprime=vdpo2(t ,y) 

W=12000; 

IXX=1760 ; IYY=9450 ; 

IZZ=10700 ; IXZ=0; IXY=0; 

IYZ=0 ; L=17.425; RH0=1.94; ZB=0; 

G=32.2; U=5 ; M=W/G; B=W; XB=0; 

ZG=0 .05 ; CD=0 . 5 ; YG=0; 

XG=1 ; YB=0; 

ND1=0.5*RH0*L‘2; 

7. DEFINE HYDRODYNAMIC COEFFICIENTS 

YPD0T=1 . 270e-04*NDl*L~2; 

YVD0T=-5 . 550e-02*NDl*L; 

YRD0T=1 . 240e-03*NDl*L~2; 
YP=3.055e-03*NDl*L; 

YV=-9 . 310e-02*NDl ; 

YR=-5 . 940e-02*NDl*L ; 

NPD0T=-3 . 370e-05*NDl*L‘3; 

NVD0T=1 . 240e-03*NDl*L~2 ; 

NRD0T=-3 .400e-03*NDl*L~3; 

NP=-8 . 405e-04*NDl*L‘2 ; 

NV=-1 ,484e-02*NDl*L; 

NR=-1 . 640e-02*NDl*L‘2 ; 

KPD0T=-1.01e-03*NDl*L~3; 

KVD0T=1.27e-04*NDl*L~2; 

KRD0T=-3 . 37e-05*NDl*L~3; 

KP=-1 . 10e-02*NDl*L‘2 ; 

KV=3.055e-03*NDl*L; 

KR=-8.41e-04*NDl*L‘2; 

D=( (M-YVDOT) *(IZZ-NRDOT)* (IXX-KPDOT) ) . . . 

- ( (M-YVDOT) * (IXZ-KRDOT) * (-IXZ-NPDOT) ) . . . 

- ( (M*XG-NVDOT) * (M*XG-YRDOT) * ( IXX-KPDOT) ) . . 
+ ( (M*XG-NVDOT) * (IXZ-KRDOT) * (-M*ZG-YPDOT) ) . 
+ ( (-M*ZG-KVDOT) * (M*XG-YRDOT) * (-IXZ-NPDOT) ) 
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- ( (-M*ZG-KVDOT) * (IZZ-NRDOT) * (-M*ZG-YPDOT) ) ; 



A1 1= ( (IZZ-NRDOT) * (IXX-KPDOT) ) - ( (IXZ-KRDOT) * (-IXZ-NPDOT) ) ; 

A12= ( (-M+XG+YRDOT) * (IXX-KPDOT) )+( (IXZ-KRDOT) * (-M*ZG-YPDOT) ) ; 

A 13= ( (M*XG-YRDOT) * (-IXZ-NPDOT) ) -( (IZZ-NRDOT) * (-M*ZG-YPDOT) ) ; 

A21= ( (-M*XG+NVDOT) * (IXX-KPDOT) )+( (-M*ZG-KVDOT) * (-IXZ-NPDOT) ) ; 

A22= ( (M-YVDOT) * (IXX-KPDOT) ) -( (-M*ZG-KVDOT) * (-M*ZG-YPDOT) ) ; 

A23= ( (-M+YVDOT) * (-IXZ-NPDOT) ) + ( (M*XG-NVDOT) * (-M*ZG-YPDOT) ) ; 

A31= ( (M*XG-NVDOT) * ( IXZ-KRDOT) ) - ( (-M*ZG-KVDOT) * (IZZ-NRDOT) ) ; 

A32= ( (-M+YVDOT) * (IXZ-KRDOT) ) + ( (-M*ZG-KVDOT) * (M*XG-YRDOT) ) ; 

A33= ( (M-YVDOT) * ( IZZ-NRDOT) ) - ( (M*XG-NVDOT) * (M*XG~YRDOT) ) ; 

% EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 

1 

F(1 , 1)=(A11*YV*U+A12*NV*U+A13*KV*U)/D; 

F ( 1 , 2 ) = (All* (YR*U-M*U) +A12* (-M*XG*U+NR*U) +A13* (M*ZG*U+KR*U) ) /D ; 

F (1 , 3)= (A11*YP*U+A12*NP*U+A13*KP*U) /D ; 

F(1 ,4)=0; 

F(2 , 1)=(A21*YV*U+A22*NV*U+A23*KV*U)/D; 

F (2 , 2 ) = (A21* (YR*U-M*U) +A22* (-M*XG*U+NR*U) +A23* (M*ZG*U+KR*U) ) /D ; 

F (2 , 3) = (A21*YP*U+A22*NP*U+A23*KP*U) /D ; 

F(2,4)=0; 

F (3 , 1 ) = (A31*YV*U+A32*NV*U+A33*KV*U) /D ; 

F (3 , 2)=(A31* (YR*U-M*U) +A32* (-M*XG*U+NR*U) +A33* (M*ZG*U+KR*U) ) /D ; 

F (3 , 3)=(A31*YP*U+A32*NP*U+A33*KP*U) /D ; 

F(3,4)=0; 

F(4, 1)=0; 

F(4,2)=0; 

F(4,3)=l; 

F(4,4)=0; 

%==================== ============================= ========== 

% CALCULATION OF THE INTEGRATION TERMS Ei 
•/,============================= ============================= 



°/. DEFINE THE LENGTH X AND HEIGHT H TERMS FOR THE INTEGRATION 
'/. ALL IN FEET. 

XL(1)=-105 .9/12; 

XL (2) =-104. 3/ 12; 

XL(3)=-99.3/12; 

XL (4) =-94. 3/12; 

XL(5)=-87.3/12; 

XL(6)=-76.8/12; 
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XL (7) =-66 . 3/12 ; 
XL(8)=-55 .8/12; 
XL(9)=72.7/12; 
XL (10) =79. 2/12; 
XL(11)=83 . 2/12 ; 
XL(12)=87.2/12; 
XL(13)=91 .2/12; 
XL(14)=95 . 2/12 ; 
XL(15)=99 . 2/12 ; 
XL(16)=101 .2/12 
XL(17)=102 . 1/ 12 
XL (18) =103. 2/ 12 



HT(1)=0; 
HT(2)=2. 28/12; 
HT(3)=8. 24/12; 
HT(4)=13. 96/12; 
HT(5)=19. 76/12; 
HT(6)=25 . 1/12; 
HT (7) =29. 36/12; 
HT(8)=31. 85/12; 
HT(9)=31. 85/12; 
HT(10)=30. 00/12 
HT(11)=27. 84/12 
HT(12)=25. 12/12 
HT(13)=21. 44/12 
HT(14)=17. 12/12 
HT(15)=12 . 0/12 ; 
HT(16)=9. 12/12; 
HT(17)=6. 72/12; 
HT(18)=0; 



for i=l:18, 

VECl(i)=HT(i)*(y(l)+XL(i)*y(2))*(abs(y(l)+XL(i)*y(2))) ; 

VEC2(i)=HT(i)*(y(l)+XL(i)*y(2))*(abs(y(l)+XL(i)*y(2)))*XL(i); 

end 



0UT=0; 
for j=l:17, 

0UT1=0 . 5* (VEC1 ( j )+VECl ( j + 1) ) * (XL( j+1) -XL( j ) ) ; 

0UT=0UT+0UT1; 

end 
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IV=CD*OUT; 



0UT=0 ; 
for j=l : 17 , 

0UT2=0 . 5* (VEC2 ( j ) +VEC2 ( j +1) ) * (XL ( j +1) -XL ( j ) ) ; 

0UT=0UT+0UT2; 

end 

IR=CD*0UT ; 



yprime=[F(l , l)*y(l)+F(l,2)*y(2)+F(l,3)*y(3)+F(l,4)*y(4)+. . . 
(All/D)*(YG*(y(3)“2+y(2)“2)-IV+(W-B)*sin(y(4)))+. . . 

(A12/D)*(IXY*y(3)~2+IYZ*y(3)*y(2)+YG*y(l)*y(2)-IR+(XG*W-XB*B)*sin(y(4)))+ 
(A13/D)*(-IXY*y(3)*y(2)-IYZ*y(2r2-M*YG*y(l)*y(3)+(YG*W-YB*B)*cos(y(4)) . . 
- (ZG*W-ZB*B) *sin(y (4) ));... 

F(2,l)*y(l)+F(2,2)*y(2)+F(2,3)*y(3)+F(2,4)*y(4)+. . . 
(A21/D)*(YG*(y(3)“2+y(2)“2)-IV+(W-B)*sin(y(4)))+. . . 

(A22/D) *(IXY*y(3) ~2+IYZ*y(3)*y(2)+YG*y (1) *y(2)-IR+(XG*W-XB*B) *sin(y(4) ) ) + 
(A23/D)*(-IXY*y(3)*y(2)-IYZ*y(2)~2-M*YG*y(l)*y(3) + (YG*W-YB*B)*cos(y(4)) . . 
- (ZG*W-ZB*B) *sin(y(4) ) ) ; . . . 

F(3,l)*y(l)+F(3,2)*y(2)+F(3,3)*y(3)+F(3,4)*y(4)+. . . 

(A31/D) *(YG*(y(3)~2+y(2)~2)-IV+(W-B) *sin(y(4) ))+. . . 

(A32/D) * (IXY*y (3) ‘2+IYZ*y (3) *y(2) +YG*y (1) *y(2)-IR+(XG*W-XB*B) *sin(y (4) ) ) + 
(A33/D)*(-IXY*y(3)*y(2)-IYZ*y(2)~2-M*YG*y(l)*y(3)+(YG*W-YB*B)*cos(y(4)) . . 
-(ZG*W-ZB*B)*sin(y(4) ));... 
l*y(3)] ; 
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HOPF . FOR 



C 
C 

C HOPF BIFURCATIONS PROGRAM 
C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

REAL*8 L , IYY , M , YPDOT , YVDOT , ND1 , YRDOT 
REAL*8 YP , YV , YR , NPDOT , NVDOT , NRDOT , NP , NV , NR 
REAL*8 GAMA , U , KVDOT , KRDOT , KPDOT , D 
REAL*8 E0,E1,E2,E3,E4,XG,ZG,KR,KP,KV,XG1,ZG1 
C 

REAL*8 Mil ,M12 ,M13 ,M14 ,M21 ,M22 ,M23 
REAL *8 M24 , M31 , M32 , M33 , M34 , M41 , M42 , M43 , M44 
REAL*8 Nil ,N12 , N13,N14 ,N21 ,N22 ,N23 ,N24 
REAL* 8 N31 , N32 , N33 , N34 , N41 , N42 , N43 , N44 
REAL* 8 LI 1 , L12 , L13 , L14 , L21 , L22 , L23 , L24 , L31 
REAL*8 L32 , L33 , L34 , L15 , L16 , L17 , L25 , L26 , L27 , L35 
REAL*8 L36 , L37 , LI A , L2A , L3A , L4A , L5A , L6A , L7A , L8A , L9A 
REAL*8 L10A,L11A , L12A, LI ,L2 ,L3, L4 ,L5 ,L6 ,L7 
REAL*8 L8 , L9 , R1 1 , R1 2 , R1 3 , R14 , R21 , R22 , R23 , R24 
REAL*8 P11,P12,P13,P21, P22 ,P23 
C 

DIMENSION F(4,4) ,T(4,4) ,TINV(4,4) ,FV1(4) , IV1 (4) ,YYY(4,4) 
DIMENSION WR(4) ,WI(4) ,TSAVE(4,4) ,TLUD(4,4) , IVLUD(4) ,SVLUD(4) 
DIMENSION ASAVE(4 ,4) ,A2(4,4) ,XL(18) ,HT(18) , ZGG(197) ,FF(4,4) 
DIMENSION VEC0C18) ,VEC1(18) ,VEC2(18) ,VEC3(18) ,VEC4(18) ,XGG(197) 
C 

INTEGER I , J,K 
C 

OPEN ( 20 , FILE= ’ HOPF . RES ’ , STATUS= > OLD > ) 

OPEN (21, FILE= ’ DATA . DAT * , STATUS= ’ OLD ’ ) 

OPEN ( 23 , FILE= * HOPF 1 . RES 9 , STATUS= * OLD ’ ) 

OPEN ( 25 , FILE= ’ H0PF2 . RES 1 , STATUS= J OLD , ) 

C 

WEIGHT=12000 . 0 
IXX =1760.0 
IYY =9450.0 
IZZ =10700.0 
IXZ =0.0 
IXY =0.0 
IYZ =0.0 
L =17.425 
RHO =1.94 
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WRITE (*,*) ’ ENTER CD’ 

READ (*,*) CD 
G =32.2 

XB =0.0 

ZB =0.0 

YG =0.0 

YB =0.0 

YDELTAR=0 . 0 
DELTAR=0 . 0 
NDELTAR=0 . 0 
NPR0P=0 . 0 
M= WEIGHT/G 
B= WEIGHT 
W=WEIGHT 
U=8 . 0 

ND1=0 . 5*RH0*L**2 
C 

C DEFINE HYDRODYNAMIC COEFFICIENTS 
C 

YPD0T=1 . 270E-04*ND1*L**2 
YVD0T=-5 . 550E-02*ND1*L 
YRD0T=1 . 240E-03*ND1*L**2 
YP=3.055E-03*ND1*L 
YV=-9 . 310E-02*ND1 
YR=-5 . 940E-02*ND1*L 
C 

NPD0T=-3 . 370E-05*ND1*L**3 
NVD0T=1.240E-03*ND1*L**2 
NRDOT=-3.400E-03*ND1*L**3 
NP=-8 . 405E-04*ND1*L**2 
NV=-1 .484E-02*ND1*L 
NR=-1 . 640E-02*ND1*L**2 
C 

KPD0T=-1 . 01E-03*ND1*L**3 
KVD0T=1 . 27E-04*ND1*L**2 
KRD0T=-3 . 37E-05*ND1*L**3 
KP=-1 . 10E-02*ND1*L**2 
KV=3 . 055E-03*ND1*L 
KR=-8 .41E-04*ND1*L**2 
C 

C DEFINE THE LENGTH X AND HEIGHT H TERMS FOR 
C THE INTEGRATION, ALL IN FEET. 



63 



c 



c 



c 



XL( 1)=-105. 9/12.0 
XL( 2)=-104. 3/12.0 
XL ( 3)=-99. 3/12.0 
XL( 4)=-94. 3/12.0 
XL( 5)=-87. 3/12.0 
XL( 6)=-76. 8/12.0 
XL ( 7)=-66. 3/12.0 
XL ( 8)=-55. 8/12.0 
XL( 9)=72. 7/12.0 
XL(10)=79. 2/12.0 
XL(11)=83. 2/12.0 
XL(12)=87. 2/12.0 
XL(13)=91. 2/12.0 
XL(14)=95. 2/12.0 
XL(15)=99. 2/12.0 
XL(16)=101. 2/12.0 
XL(17)=102. 1/12.0 
XL(18)=103. 2/12.0 

HT( 1)= 0.000 
HT( 2)= 2.28/12.0 
HT( 3)= 8.24/12.0 
HT( 4)= 13.96/12.0 
HT ( 5)= 19.76/12.0 
HT ( 6)= 25.1/12.0 
HT ( 7)= 29.36/12.0 
HT ( 8) = 31.85/12.0 
HT ( 9)= 31.85/12.0 
HT(10)= 30.00/12.0 
HT ( 1 1 ) = 27.84/12.0 
HT(12)= 25.12/12.0 
HT(13)= 21.44/12.0 
HT(14)= 17.12/12.0 
HT(15)= 12.0/12.0 
HT(16)= 9.12/12.0 
HT(17)= 6.72/12.0 
HT(18)= 0.00 

DO 104 K = 1,18 
VECO(K)=HT (K) 
VEC1(K)=XL(K)*HT(K) 
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VEC2 (K) =XL (K) *XL (K) *HT (K) 

VEC3 (K) =XL (K) *XL (K) *XL (K) *HT (K) 

VEC4 (K) =XL (K) *XL (K) *XL (K) *XL (K) *HT (K) 

104 CONTINUE 

CALL TRAP ( 18 , VECO , XL , EO ) 

CALL TRAP ( 18 , VEC1 , XL , El ) 

CALL TRAP(18,VEC2,XL,E2) 

CALL TRAP ( 18 , VEC3 , XL , E3) 

CALL TRAP (18, VEC4 , XL , E4) 

C 

WRITE (*,*) ’ ENTER GAMA’ 

READ (*,*) GAMA 

C READ THE CRITICAL VALUES FOR XG AND ZG FROM FILE DATA . DAT 

C 

XGG(1)=0.0 
ZGG(1)=0. 016358083 
DO 1 1=2,197 

READ (21 , *)XG ,ZG 

XGG(I)=XG 

ZGG(I)=ZG 

C DETERMINE [F] COEFFICIENTS 

C 

D= ( (M-YVDOT) * ( I ZZ-NRDOT) * ( IXX-KPDOT) ) 

& - ( (M-YVDOT) * (IXZ-KRDOT) * (-IXZ-NPDOT) ) 

& - ( (M*XG-NVDOT) * (M*XG-YRDOT) * (IXX-KPDOT) ) 

& + ( (M*XG-NVDOT) * ( IXZ-KRDOT) * ( -M*ZG-YPDOT) ) 

& + ( (-M*ZG-KVDOT) * (M*XG-YRDOT) * ( -IXZ-NPDOT) ) 

& - ( (-M*ZG-KVDOT) * (IZZ-NRDOT) * (-M*ZG-YPDOT) ) 

C 

All=( (IZZ-NRDOT) * (IXX-KPDOT) )-( (IXZ-KRDOT) *(-IXZ-NPDOT) ) 
A12=( (-M*XG+YRDOT) * (IXX-KPDOT) )+( (IXZ-KRDOT) * (-M*ZG-YPDOT) ) 
A13= ( (M*XG-YRDOT) * (-IXZ-NPDOT) ) - ( ( IZZ-NRDOT) * (-M*ZG-YPDOT) ) 
A21= ( (-M*XG+NVDOT) * (IXX-KPDOT) ) + ( (-M*ZG-KVDOT) * (-IXZ-NPDOT) ) 
A22=( (M-YVDOT) * (IXX-KPDOT) ) - ( (-M*ZG-KVDOT) * (-M*ZG-YPDOT) ) 
A23= ( (-M+YVDOT) * (-IXZ-NPDOT) ) + ( (M*XG-NVDOT) * (-M*ZG-YPDOT) ) 
A31=( (M*XG-NVDOT) * (IXZ-KRDOT) ) - ( (-M*ZG-KVDOT) * (IZZ-NRDOT) ) 
A32=( (-M+YVDOT) * (IXZ-KRDOT) )+ ( (-M*ZG-KVDOT) * (M*XG-YRDOT) ) 
A33= ( (M-YVDOT) * (IZZ-NRDOT) ) - ( (M*XG-NVDOT) * (M*XG-YRDOT) ) 

C 

F(1,1)=(A11*YV*U+A12*NV*U+A13*KV*U)/D 
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F(1,2)=(A11*(YR*U-M*U)+A12*(-M*XG*U+NR*U)+A13*(M*ZG*U+KR*U))/D 
F(1,3)=(A11*YP*U+A12*NP*U+A13*KP*U)/D 
F(1,4)=(A11*(W-B)+A12*(XG*W-XB*B)+A13*(-ZG*W+ZB*B))/D 
F (2 , 1) = (A2 1*YV*U+A22*NV*U+A23*KV*U) /D 

F (2 , 2) = (A21* (YR*U-M*U)+A22*(-M*XG*U+NR*U) +A23* (M*ZG*U+KR*U) ) /D 
F (2 , 3) = (A21*YP*U+A22*NP*U+A23*KP*U) /D 
F(2,4)=(A21*(W-B)+A22*(XG*W-XB*B)+A23*(-ZG*W+ZB*B) )/D 
F(3 , 1)=(A31*YV*U+A32*NV*U+A33*KV*U) /D 

F(3,2)=(A31* (YR*U-M*U)+A32* (-M*XG*U+NR*U) +A33* (M*ZG*U+KR*U) ) /D 
F(3 , 3)=(A31*YP*U+A32*NP*U+A33*KP*U) /D 
F(3,4)=(A31*(W-B)+A32*(XG*W-XB*B)+A33*(-ZG*W+ZB*B))/D 
F (4,1) =0.0 
F(4,2)=0.0 
F(4,3)=l .0 
F (4,4)=0 . 0 
C 

C EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 

C 

DO 11 K=1 , 4 
DO 12 J=1 , 4 

ASAVECK, J)=F(K, J) 

12 CONTINUE 
11 CONTINUE 

CALL RG (4,4,F, WR, WI , 1 , YYY , IV1 ,FV1 , IERR) 

WRITE(23 , 1007)WR(1) ,WR(2) ,WR(3) ,WR(4) 

CALL DSOMEG ( IEV , WR , WI , OMEGA , CHECK) 

OMEGAO=OMEGA 
DO 5 J=1 ,4 

T(J, 1)= YYY(J.IEV) 

T( J , 2) =-YYY( J , IEV+1) 

5 CONTINUE 

IF (IEV. EQ. 1.0) GO TO 13 
IF (IEV. EQ. 2.0) GO TO 18 
IF (IEV. EQ. 3.0) GO TO 14 
STOP 3004 
14 DO 6 J=1 ,4 

T( J ,3) =YYY (J , 1) 

T( J ,4) =YYY ( J , 2) 

6 CONTINUE 
GO TO 17 

18 DO 19 J=1 , 4 

T( J , 3) =YYY ( J , 1) 
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T( J , 4)=YYY ( J ,4) 

19 CONTINUE 
GO TO 17 

13 DO 16 J=l,4 

T ( J , 3)=YYY( J ,3) 

T ( J ,4)=YYY ( J ,4) 

16 CONTINUE 

17 CONTINUE 
C 

C NORMALIZATION OF THE CRITICAL EIGENVECTOR 

C 

CALL NORMAL (T) 

C 

C INVERT TRANSFORMATION MATRIX 

C 

DO 2 K=1 ,4 
DO 3 J=1 ,4 
TINV(K, J)=0.0 
TSAVECK, J)=T(K,J) 

3 CONTINUE 
2 CONTINUE 

CALL DLUD (4,4, TSAVE , 4 , TLUD , IVLUD) 

DO 4 J=l,4 

IF (IVLUD ( J) .EQ.O) STOP 3003 

4 CONTINUE 

CALL DILUC4, 4, TLUD, IVLUD, SVLUD) 

DO 8 K=1 ,4 
DO 9 J=1 ,4 

TINVCK, J)=TLUD(K, J) 

9 CONTINUE 

8 CONTINUE 
C 

C CHECK Inv(T) *A*T 

C 

CALL MULT (TINV , ASAVE , T , A2) 

C 

P1=A2(3,3) 

P2=A2(4,4) 

PEIG1=P1 

PEIG2=P2 

WRITE(25 , 1008)A2(1 , 1) , A2 (2, 2) , PI , P2 
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c 

c 



c 

c 

c 



c 

c 

c 



DEFINITION OF Nij 

N11=TINV(1 , 1) 

N12=TINV (1,2) 

N13=TINV(1,3) 

N14=TINV(1 ,4) 

N21=TINV (2 ,1) 

N22=TINV(2,2) 

N23=TINV (2,3) 

N24=TINV (2,4) 

N31=TINV(3,1) 

N32=TINV (3,2) 

N33=TINV(3,3) 

N34=TINV (3,4) 

N41=TINV (4,1) 

N42=TINV (4 ,2) 

N43=TINV(4,3) 

N44=TINV (4,4) 

DEFINITION OF Mij 

M11=T(1 , 1) 

M12=T(1,2) 

M13=T(1 ,3) 

M14=T(1 ,4) 

M21=T(2 ,1) 

M22=T(2,2) 

M23=T(2,3) 

M24=T(2 ,4) 

M31=T (3,1) 

M32=T(3,2) 

M33=T (3 ,3) 

M34=T(3,4) 

M41=T(4, 1) 

M42=T(4,2) 

M43=T (4 , 3) 

M44=T (4,4) 

DEFINITION OF Lij 

L1=YG*((M31**2)+(M21**2)) 

L2=IXY* (M31**2)+IYZ*M31*M21+YG*M21*M11 
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L3=IXY*M31*M21+IYZ*(M21**2)+M*YG*M11*M31+YG*W*(M41**2) 

& -YB*B* (M41**2) 

L4=YG* (2 . 0*M31*M32+2 . 0*M21*M22) 

L5=2.0*IXY*M31*M32+IYZ*(M31*M22+M32*M21)+YG*(M21*M12+M22*M11) 
L6=IXY*(M31*M22+M32*M21)+2.0*IYZ*M21*M22+M*YG* (M11*M32+M12*M31) 

& +2 . 0* (YG*W-YB*B) *M41*M42 

L7=YG* ((M32**2)+(M22**2)) 

L8=IXY* (M32**2)+IYZ*M32*M22+YG*M22*M12 

L9=IXY*M32*M22+IYZ* (M22**2)+M*YG*M12*M32+(YG*W-YB*B)* (M42**2) 

C 

C 

L15=(A11/D)*L1+(A12/D)*L2-(A13/D)*L3 
L16=(A11/D)*L4+(A12/D)*L5-(A13/D)*L6 
L17=(A11/D)*L7+(A12/D)*L8-(A13/D)*L9 
L25=( A21/D) *L1+(A22/D)*L2- (A23/D) *L3 
L26=(A21/D)*L4+(A22/D)*L5-(A23/D)*L6 
L27=(A21/D)*L7+(A22/D)*L8-(A23/D) *L9 
L35=(A31/D) *L1+(A32/D)*L2-(A33/D)*L3 
L36= ( A31/D) *L4+(A32/D)*L5- (A33/D) *L6 
L37=(A31/D)*L7+(A32/D)*L8- (A33/D) *L9 
C 

C=CD/(6.0*GAMA) 

L1A=C*(E0*(M11**3)+3.0*E1*(M11**2)*M21+3.0*E2*M11*(M21**2) 

& +E3*(M21**3))+(1.0/6.0)*(W-B)*(M41**3) 

L2A=C* (El* (Mll**3)+3 . 0*E2* (Mll**2) *M21+3 . 0*E3*M11* (M21**2) 

& +E4* (M21**3) ) + ( 1 . 0/6 . 0) * (XG*W-XB*B) * (M41**3) 

L3A= ( 1 . 0/6 . 0) * (ZG*W-ZB*B) * (M41**3) 

L4A=C*(3.0*E0* (M 1 1**2) *M12+3.0*E1*( (Ml 1**2) *M22+2 .0*M11*M12*M21) 
& +3 . 0*E2* (M12* (M21**2)+2 . 0*M21*M22*Mll)+3 .0*E3* (M21**2) *M22) 

& +(1 .0/6.0) *(W-B) *3.0* (M41**2)*M42 

L5A=C* (3 . 0*E1* (Mll**2)*M12+3 . 0*E2* ( (Mll**2) *M22+2 . 0*M1 1*M12*M21) 
& +3.0*E3*(M12*(M21**2)+2.0*M21*M22*M11)+3.0*E4*(M21**2)*M22) 

& + ( 1 . 0/6 . 0) * (XG*W-XB*B) *3 . 0* (M41**2) *M42 

L6A=(1.0/6.0)*(ZG*W-ZB*B)*3.0*(M41**2)*M42 

L7A=C*(3.0*E0*M11*(M12**2)+3.0*E1*((M12**2)*M21+2.0*M11*M12*M22) 
& +3 . 0*E2* ( (M22**2) *M1 1+2 . 0*M21*M22*M12)+3 . 0*E3*M21* (M22**2) ) 

& +(1.0/6.0)*(W-B)*3.0*M41*(M42**2) 

L8A=C* (3 .0*E1*M11* (M12**2)+3 . 0*E2* ( (M12**2) *M21+2 . 0*M11*M12*M22) 
& +3.0*E3*((M22**2)*M11+2.0*M21*M22*M12)+3.0*E4*M21*(M22**2)) 

& + ( 1 . 0/6 . 0) * (XG*W-XB*B) *3 . 0*M41* (M42**2) 

L9A= ( 1 . 0/6 . 0) * (ZG*W-ZB*B) *3 . 0*M41* (M42**2) 

L10A=C* (E0* (M12**3)+3 . 0*E1* (Mll**2) *M21+3 . 0*E2*M12* (M22**2) 
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& +E3*(M22**3))+(1.0/6.0)*(W-B)*(M42**3) 

L11A=C* (El* (Ml 2**3) +3 . 0*E2* (Ml 1**2) *M21+3 . 0*E3*M12* (M22**2) 
& +E4* (M22**3) )+(l .0/6.0)* (XG*W-XB*B) * (M42**3) 

L12A=(1.0/6.0)*(ZG*W-ZB*B)*(M42**3) 

C 

Lll= (-A11/D) *LlA+(-A12/D) *L2A+(A13/D) *L3A 
L12= (-A1 1/D) *L4A+(-A12/D) *L5A+ (A13/D) *L6A 
LI 3= (- All/D) *L7A+(-A12/D) *L8A+(A13/D) *L9A 
L14= (- All/D) *L10A+(-A12/D) *L11A+(A13/D) *L12A 
C 

L21=(-A21/D) *LlA+(-A22/D) *L2A+(A23/D)*L3A 
L22=(-A21/D)*L4A+(-A22/D)*L5A+(A23/D)*L6A 
L23= (-A21/D) *L7A+(-A22/D) *L8A+(A23/D) *L9A 
L24=(-A21/D)*L10A+(-A22/D)*LllA+(A23/D)*L12A 
C 

L31=(-A31/D) *LlA+(-A32/D) *L2A+(A33/D) *L3A 
L32=(-A31/D) *L4A+(-A32/D) *L5A+(A33/D) *L6A 
L33=(-A31/D)*L7A+(-A32/D)*L8A+(A33/D)*L9A 
L34= (-A31/D) *L10A+(-A32/D) *L11A+(A33/D) *L12A 
C 

R11=(N11*L11)+(N12*L21)+(N13*L31) 

R12= (N11*L12)+(N12*L22) +(N13*L32) 
R13=(N11*L13)+(N12*L23)+(N13*L33) 

R14= (N11*L14)+(N12*L24)+(N13*L34) 
R21=(N21*L11)+(N22*L21)+(N23*L31) 

R22= (N21*L12)+(N22*L22)+(N23*L32) 

R23=(N21*L13)+(N22*L23) +(N23*L33) 

R24=(N21*L14) + (N22*L24) + (N23*L34) 

C 

P11=(N11*L15)+(N12*L25)+(N13*L35) 

P12=(N11*L16)+(N12*L26)+(N13*L36) 

P13=(N11*L17)+(N12*L27)+(N13*L37) 

P21=(N21*L15)+(N22*L25) +(N23*L35) 
P22=(N21*L16)+(N22*L26)+(N23*L36) 

P23= (N21*L17) + (N22*L27 ) + (N23*L37 ) 

C 

C EVALUATE DALPHA AND DOMEGA 

C 

ZGR =ZGG(I) 

ZGL =ZGG(I-1) 

ZG1 =ZGR 
XG1 =XGG(I) 
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c 

CALL FMATRIX(ZG1,XG1,FF) 

C 

CALL RG(4,4,FF,WR,WI ,0, YYY, IV1 ,FV1 , IERR) 

CALL DSTABL(DEOS,WR,WI , FREQ) 

ALPHR=DE0S 

0MEGR=FREQ 

C 

ZG1 =ZGL 
XG1 =XGG(I) 

C 

CALL FMATRIX (ZG1 , XG1 , FF) 

C 

CALL RG(4,4, FF , WR, WI , 0 , YYY, IV1 , FV1 , IERR) 

CALL DSTABL (DEOS , WR , WI , FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 

C 

DALPHA= (ALPHR-ALPHL) / (ZGR-ZGL) 

DOMEGA= (OMEGR-OMEGL) / (ZGR-ZGL) 

C 

C EVALUATION OF HOPF BIFURCATION COEFFICIENTS 

C 

C0EF1= (1 . 0/8 . 0) * (3 . 0*R1 1+R13+R22+3 . 0*R24) 

C0EF2= (1 . 0/8 . 0) * (3 . 0*Rll+R23-R12-3 . 0*R14) 

WRITE (20 , 2001 ) XG , ZG , C0EF1 , DALPHA , OMEGAO , PEIG1 , PEIG2 
1 CONTINUE 
C 

STOP 

2001 FORMAT (7E14.5) 

1007 FORMAT (4E14.5) 

1008 FORMAT (4E14.5) 

END 
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